#include <set>
#include "misclib.h"
#include "ssystem.h"
#include "model.h"
#include "datatypes.h"
#ifndef SERIALQGCM
#include <omp.h> // include openMPI
#endif

#ifdef DEBUG
#define dbgprintf(arg,...) mo.printfmsg(arg, ## __VA_ARGS__)
#else
#define dbgprintf(arg,...) // just eliminate the code if we are compiling an optimized version
#endif

qgcomplex null_complex(0,0);
/* ------------------------------------------------------------------------------
   lopinfo class: used to store information about the linear operators
   also defined: lopinfo::element (an element of a  linear operator)
   -------------------------------------------------------------------------------- */
lopinfo::element::element(){
	field_indices = rank2(-1,-1);
	angular_mode = -1;
	parser = NULL;
}
// - - - - - - - - - - - -
lopinfo::element::element(const lopinfo::element & ol){
	field_indices = ol.field_indices;
	angular_mode = ol.angular_mode;
	if(ol.parser != NULL ){
		parser = new ParserX(*ol.parser);
	}else{
		parser = NULL;
	}
}
// - - - - - - - - - - - -
lopinfo::element & lopinfo::element::operator=(const lopinfo::element & ol){
	field_indices = ol.field_indices;
	angular_mode = ol.angular_mode;
	if(ol.parser != NULL ){
		parser = new ParserX(*ol.parser);
	}else{
		parser = NULL;
	}
	return *this;
}
// - - - - - - - - - - - -
lopinfo::element::~element(){
	delete this->parser;
}
// - - - - - - - - - - - -
void lopinfo::element::SetExpr(const token_list & expr){
	if(parser == NULL) parser = new ParserX;
	parser->SetExpr(expr);
	try {
		parser->Eval();
	}catch(const ParserError & e){
		fprintf(stderr, "Parser evaluation failed with message: %s\n", e.what());
		fatal_stop("Parser evaluation failed.");
	}
}
// - - - - - - - - - - - -
REAL lopinfo::element::imag() const{
	rassert( parser != NULL );
	return static_cast<float>( parser->Eval().GetComplex().imag() );
};
// - - - - - - - - - - - -
REAL lopinfo::element::real() const{
	rassert( parser != NULL );
	return static_cast<float>( parser->Eval().GetComplex().real() );
};
// - - - - - - - - - - - -
void lopinfo::element::assign(const INDEX & angm, const rank2 & fpos, const token_list & expr, const ParserX & Pinit){
	angular_mode = angm;
	field_indices = fpos;
	if(parser == NULL) parser = new ParserX(Pinit);
	parser->SetExpr(expr);
};
// - - - - - - - - - - - -
lopinfo::lopinfo(){}
lopinfo::~lopinfo(){}

// ----------------------------------------------------------------------------- //
// netcdf wrapper/auxiliary functions (written locally here)
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void nc_assert_att_exists(NcAtt * att, NcVar * var, const char * att_name){
	if(att) return;
	fprintf(stderr, "Error: Could not retrieve attribute %s from %s.\n", att_name, var->name());
	NcAtt * atype = var->get_att("type");
	if(atype){
		char * type_name =  atype->as_string(0);
		fprintf(stderr, "Error: Unable to read %s data.\n", type_name);
		delete [] type_name;
	}
	fprintf(stderr, "       Is the model state file corrupt?\n");
	throw false;
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void assert_no_nc_error( NcError & error_status, NcVar * var, const char * att_name ){
	int err = error_status.get_err();
	if( err != NC_NOERR ){
		fprintf(stderr, "NETCDF Error: %s.\n", error_status.get_errmsg());
		// additional information
		fprintf(stderr, "Error reading %s.\n", var->name());
		NcAtt * atype = var->get_att("type");
		if(atype){
			char * type_name =  atype->as_string(0);
			fprintf(stderr, "Error: Unable to read %s data.\n", type_name);
			delete [] type_name;
		}
		fprintf(stderr, "       Is the model state file corrupt?\n");
		throw false;
	}
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
int get_integer_att(NcVar * var, const char * name){
	NcError error_status( NcError::silent_nonfatal );
	NcAtt * att = var->get_att(name);
	nc_assert_att_exists( att , var, name ); 
	int rvalue = att->as_int(0);
	assert_no_nc_error( error_status, var, name );
	return rvalue;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
string get_string_att(NcVar * var, const char * name){
	NcError error_status( NcError::silent_nonfatal );
	NcAtt * att = var->get_att(name);
	nc_assert_att_exists( att, var, name );
	char * tmp_rvalue = att->as_string(0);
	assert_no_nc_error( error_status, var, name );
	string rvalue(tmp_rvalue); delete [] tmp_rvalue;
	return rvalue;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void check_dissip_ops(NcFile * ncfile, const vector< pair<INDEX, ParserX> > & dissipation_ops){
	NcError error_status( NcError::silent_nonfatal );
	// load the dissipation operators
	int iter = 0;
	vector<bool> dissip_op_check_list;
	if( dissipation_ops.size() > 0 ) dissip_op_check_list.resize( dissipation_ops.size(), false );
	while(true){
		const string nc_var_name = string("dissip_op_") + convertint2string(iter);
		NcVar * var = ncfile->get_var(nc_var_name.c_str());
		if( var == NULL ) break;
		const int findex = get_integer_att(var, "field_index");
		const string op_expr = get_string_att(var, "expr");
		// we have a dissipation operator, they don't always have names, so we do our best to match up
		// existing with previous data, there is no requirement that all operators exist nor that
		// previously defined operators exist ...
		bool found_match = false;
		for(size_t u = 0; u < dissipation_ops.size(); u++){
			if( dissipation_ops[u].first == static_cast<INDEX>(findex) && 
					dissipation_ops[u].second.GetExprAsString() == op_expr ){   
				// we probably have a match
				found_match = true;
				dissip_op_check_list[u] = true;
				break;
			}
		}
		if( ! found_match ){ // print a warning
      fprintf(stderr, "Warning: Dissipation operator #%i in saved model state has no matching definition\n", iter);
      fprintf(stderr, "         in current model instantiation.\n");
      fprintf(stderr, "         Operator information:\n");
      fprintf(stderr, "          Present in the model's saved state.\n");
      fprintf(stderr, "           field index: '%i'\n", findex+1);
      fprintf(stderr, "           expression:  '%s'\n", op_expr.c_str());
    }
    iter++;
  }
  // now check check_list to make sure all instantiated operators were matched by the model's saved state
  for(size_t u = 0; u < dissip_op_check_list.size(); u++){
    if(!dissip_op_check_list[u]){
      fprintf(stderr, "Warning: Dissipation operator #%i defined in the current instantiation of the model\n"
          , static_cast<int>(u));
      fprintf(stderr, "         has no matching definition in the saved model state.\n");
      fprintf(stderr, "         Operator information:\n"); 
      fprintf(stderr, "          Present in the model's saved state.\n");
      fprintf(stderr, "           field index: '%i'\n", static_cast<int>(dissipation_ops[u].first)+1);
      fprintf(stderr, "           expression:  '%s'\n", dissipation_ops[u].second.GetExprAsString().c_str());
    }
  } 
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void check_inhomogeneous_ops(NcFile * ncfile, const vector< pair<rank3, ParserX> > & inhomogeneous_ops){
  NcError error_status( NcError::silent_nonfatal );
  // load the inhomogeneous operators
  int iter = 0;
  vector<bool> ihomo_op_check_list;
  if( inhomogeneous_ops.size() > 0 ) ihomo_op_check_list.resize( inhomogeneous_ops.size(), false );
  while(true){
    const string nc_var_name = string("inhomogeneous_op_") + convertint2string(iter);
    NcVar * var = ncfile->get_var(nc_var_name.c_str());
    if( var == NULL ) break;
    const int nkindex = get_integer_att(var, "nkindex");
    const int Iindex = get_integer_att(var, "Iindex");
    const int Jindex = get_integer_att(var, "Jindex");
    const string op_expr = get_string_att(var, "expr");
    const rank3 op_pos(static_cast<INDEX>(nkindex), static_cast<INDEX>(Iindex), static_cast<INDEX>(Jindex));
    // we have a inhomogeneous operator, they don't have names, so we do our best to match up
    // existing with previous data, there is no requirement that all operators exist nor that
    // previously defined operators exist ...
    bool found_match = false;
    for(size_t u = 0; u < inhomogeneous_ops.size(); u++){
      const rank3 & r3 = inhomogeneous_ops[u].first; 
      if( ( r3 == op_pos ) && 
          inhomogeneous_ops[u].second.GetExprAsString() == op_expr ){   
        // we probably have a match
        found_match = true;
        ihomo_op_check_list[u] = true;
        break;
      }
    }
    if( ! found_match ){ // print a warning
      fprintf(stderr, "Warning: Inhomogeneous operator #%i in saved model state has no matching definition\n", iter);
      fprintf(stderr, "         in current model instantiation.\n");
      fprintf(stderr, "         Operator information:\n");
      fprintf(stderr, "          Present in the model's saved state.\n");
      fprintf(stderr, "           nki=%i\ti=%i\tj=%i", nkindex, Iindex+1, Jindex+1);
      fprintf(stderr, "           expression:  '%s'\n", op_expr.c_str());
    }
    iter++;
  }
  // now check check_list to make sure all instantiated operators were matched by the model's saved state
  for(size_t u = 0; u < ihomo_op_check_list.size(); u++){
    if(!ihomo_op_check_list[u]){
      const rank3 & r3 = inhomogeneous_ops[u].first;
      fprintf(stderr, "Warning: Inhomogeneous operator #%i defined in the current instantiation of the model\n", 
          static_cast<int>(u));
      fprintf(stderr, "         has no matching definition in the saved model state.\n");
      fprintf(stderr, "         Operator information:\n"); 
      fprintf(stderr, "          Present in the model's saved state.\n");
      fprintf(stderr, "           nki=%i\ti=%i\tj=%i", static_cast<int>(r3[0]), static_cast<int>(r3[1])+1,static_cast<int>(r3[2])+1);
      fprintf(stderr, "           expression:  '%s'\n", inhomogeneous_ops[u].second.GetExprAsString().c_str());
    }
  } 
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void check_linear_ops(NcFile * ncfile, const vector<lopinfo::element> & linear_ops){
  NcError error_status( NcError::silent_nonfatal );
  // load the linear operators
  int iter = 0;
  vector<bool> linear_op_check_list;
  if( linear_ops.size() > 0 ) linear_op_check_list.resize( linear_ops.size(), false );
  while(true){
    const string nc_var_name = string("linear_op_") + convertint2string(iter);
    NcVar * var = ncfile->get_var(nc_var_name.c_str());
    if( var == NULL ) break;
    const int nkindex = get_integer_att(var, "nkindex");
    const int Iindex = get_integer_att(var, "Iindex");
    const int Jindex = get_integer_att(var, "Jindex");
    const string op_expr = get_string_att(var, "expr");
    const rank3 op_pos(static_cast<INDEX>(nkindex), static_cast<INDEX>(Iindex), static_cast<INDEX>(Jindex));
    // we have a linear operator, they don't have names, so we do our best to match up
    // existing with previous data, there is no requirement that all operators exist nor that
    // previously defined operators exist ...
    bool found_match = false;
    for(size_t u = 0; u < linear_ops.size(); u++){
      const rank3 r3(linear_ops[u].angular_mode, linear_ops[u].field_indices[0], linear_ops[u].field_indices[1]);
      if( ( r3 == op_pos ) && 
          ( linear_ops[u].parser->GetExprAsString() == op_expr )  ){
        // we probably have a match
        found_match = true;
        linear_op_check_list[u] = true;
        break;
      }
    }
    if( ! found_match ){ // print a warning
      fprintf(stderr, "Warning: Linear operator #%i in saved model state has no matching definition\n", iter);
      fprintf(stderr, "         in current model instantiation.\n");
      fprintf(stderr, "         Operator information:\n");
      fprintf(stderr, "          Present in the model's saved state.\n");
      fprintf(stderr, "           nki=%i\ti=%i\tj=%i\n", nkindex, Iindex+1, Jindex+1);
      fprintf(stderr, "           expression:  '%s'\n", op_expr.c_str());
    }
    iter++;
  }
  // now check check_list to make sure all instantiated operators were matched by the model's saved state
  for(size_t u = 0; u < linear_op_check_list.size(); u++){
    if(!linear_op_check_list[u]){
      const rank3 r3(linear_ops[u].angular_mode, linear_ops[u].field_indices[0], linear_ops[u].field_indices[0]);
      fprintf(stderr, "Warning: Linear operator #%i defined in the current instantiation of the model\n", static_cast<int>(u));
      fprintf(stderr, "         has no matching definition in the saved model state.\n");
      fprintf(stderr, "         Operator information:\n"); 
      fprintf(stderr, "          Present in the model's saved state.\n");
      fprintf(stderr, "           nki=%i\ti=%i\tj=%i\n", static_cast<int>(r3[0]), static_cast<int>(r3[1])+1,static_cast<int>(r3[2])+1);
      fprintf(stderr, "           expression:  '%s'\n", linear_ops[u].parser->GetExprAsString().c_str());
    }
  } 
};

/* ------------------------------------------------------------------------------------
   definition of the initialdata class, an intermediate data storage 
   class used to generate the initial conditions for the S-system.
   -------------------------------------------------------------------------------------- */
initialdata::initialdata(){
  index_i=Value(static_cast<int_type>(-1));
  index_j=Value(static_cast<int_type>(-1));
  index_k=Value(static_cast<int_type>(-1));
  index_mode=Value(static_cast<int_type>(-1));
  kval = Value( static_cast<float_type>(0.0) );
}
// -------------------------------------------------------------------------------
initialdata::~initialdata(){
  this->rules.clear();
}
// -------------------------------------------------------------------------------
initialdata * spreprocdata::getdata(){
  return this->initial_conditions;
}
// -------------------------------------------------------------------------------
initialdata::rule::rule(){
  // parser elements are defined by the namelist structure later
};
// -------------------------------------------------------------------------------
initialdata::rule::~rule(){
  //  reset unions
  namelist.clear();
  line_number = 0;

};
/* ------------------------------------------------------------------------------------
   Container class for all the initial_conditions rules:
   -------------------------------------------------------------------------------------- */
spreprocdata::spreprocdata(){
  this->initial_conditions = new initialdata;
};
// -------------------------------------------------------------------------------
void spreprocdata::clear(){
  this->OPdata.clear();
  delete this->initial_conditions;
  this->initial_conditions = NULL;
}
// -------------------------------------------------------------------------------
spreprocdata::~spreprocdata(){
  delete this->initial_conditions;
};
// -------------------------------------------------------------------------------
const vector<linearOPtranslationdata> & spreprocdata::getOPdata() const{
  return this->OPdata;
};
/* ------------------------------------------------------------------------------------
   Definition of local S-system class (internal stuff not needed at global scope) 
   -------------------------------------------------------------------------------------- */
class ssystem::clocal {
  public:
    // base data
    model_class & mo;
    ssystem & ssys;
    const vector<int> & nkdomain;
    const vector<REAL> & kdomain_centered;
    const vector<REAL> & dkdomain;
    const vector<REAL> & kdomain_centered_squared;
    unsigned int theta_integrated_lop_size;

    const int num_nks;
    const int num_ks;
    const int numfields;

    const bool useregularization; 
    const bool nonlinear;
    bool abort; // abort current integration and rewind (1) if flag is set
    bool initialized;
    bool stepequaltimeonly;
    bool check_realizability;
    REAL fourier_transform_coefficient;
    unsigned int tindex;
    INDEX * kminimum; // holds snapped index of alpha * kalpha for all the different ks

    vector<int> * tprimeintegrationlist;

#ifdef DEBUG // store widths (i.e. TIMESTEP associated with going from tindex to tindex+1
    vector<REAL> timestepwidths;
#else
    REAL * timestepwidths;
#endif

    twotimevar * phi_data;
    twotimevar * G_data;

    twotimevar * R_data;
    twotimevar * R_data_regularized;
    twotimevar * S_data;
    twotimevar * S_data_regularized;


    twotimevar * F_phi_data;
    twotimevar * F_phi_data_regularized;
    twotimevar * F_G_data;

    struct {
      FILE * L_out;
      FILE * R_out;
      FILE * S_out;
      FILE * PHI_out;
      FILE * FPHI_out;
      FILE * G_out;
      FILE * FG_out;
    } dump_file_handles;


  private: 
    clocal(const clocal &);
    clocal & operator=(const clocal &); 
  public:
    clocal(model_class & mdl, ssystem & S, const int nf, const vector<int> & nkd, const vector<REAL> & kc, 
        const vector<REAL> & dk, const vector<REAL> & k2,  bool nl, bool regularized) :
      mo(mdl), ssys(S), nkdomain(nkd), kdomain_centered(kc), dkdomain(dk), kdomain_centered_squared(k2), 
      num_nks(nkd.size()), num_ks(kc.size()), numfields(nf), nonlinear(nl), 
      useregularization(regularized)
  {
    this->stepequaltimeonly = false;
    this->initialized = false;
    this->tindex = 0;
    this->abort = false;
    this->tprimeintegrationlist = NULL;
#ifndef DEBUG
    this->timestepwidths = NULL;
#endif
    phi_data = NULL;
    G_data = NULL;
    R_data = NULL;
    S_data = NULL;
    S_data_regularized = NULL;
    F_phi_data = NULL;
    F_phi_data_regularized = NULL;
    F_G_data = NULL;
    if(useregularization){
      const REAL & alpha = mo.get_global_param_by_name(gdb_str_KCUTOFFMULTIPLIER).getfValue();
      this->kminimum = new INDEX[ kdomain_centered.size() ];
      for(size_t i = 0; i < kdomain_centered.size(); i++){
        REAL kmin = kdomain_centered[i] * alpha;
        int x = 0;
        if( kmin < kdomain_centered[0] ){
          kminimum[i] = -1;
        }else{
          while( kdomain_centered[x] < kmin && x < static_cast<int>(kdomain_centered.size()) ) x++;
          kminimum[i] = x;  
        }
      }
    }
  }
    // .............................................................. 
    // class ssystem_registers
    // ..............................................................
    struct ssystem_registers {
      // this structure contains variables that are accessed during the iterative evaluation
      // of the linear operators
      private:
        Value k_reg;
        Value i_reg; // field index i
        Value j_reg; // field index j
        Value m_reg; // angular mode index m
        ssystem_registers(const ssystem_registers &);
      public:
        ssystem_registers() {};
        ~ssystem_registers() {};
        void set_reg(const REAL & k)
        {
          this->k_reg = static_cast<float_type>(k);
        }; // theta_index = 0.01*theta/(2*PI)
        void set_ang(const int & m)
        {
          this->m_reg = static_cast<int_type>(m);
        }
        void set_fields(const INDEX & i, const INDEX & j)
        {
          this->i_reg = static_cast<int_type>(i);
          this->j_reg = static_cast<int_type>(j);
        }
        Variable get_k() { return Variable(&k_reg); };
        Variable get_m() { return Variable(&m_reg); };
        Variable get_i() { return Variable(&i_reg); };
        Variable get_j() { return Variable(&j_reg); };
    } reg;
    // ..............................................................
    // we define the lopinfo struct to store the information 
    // needed to evaluate each linear op (i.e. Lij(k, theta) )
    // ..............................................................
    lopinfo linear_ops; 
    vector< pair< INDEX, ParserX > > dissipation_ops; // diagonal in field index
    vector< pair< rank3, ParserX > > inhomogeneous_ops; 
    // explicit function of k, list iterates over mode, field index i, field index j 

    // ...............................................................
    // Phi storage for a 1 set of t, t' (all m, k, i, j)
    // ...............................................................
    regset phi_registers;
    regset G_registers;

    regset R_registers; // for nonlinear integration
    regset R_registers_regularized;
    regset S_registers; // for nonlinear integration
    regset S_registers_regularized;

    regset F_phi_registers;       // RHS of d(phi)/dt unequal time including linear part
    regset F_phi_registers_regularized;

    regset F_G_registers;       // RHS of d(G)/dt   unequal time

    regset lop_register;
    regset lop_non_inverted; // before multiplying by QI
    regset dissip_register;
    regset inhomo_forcing_register; // inhomogeneous forcing term
    // . ...................................................................................
    // extended position storage class 
    // . ...................................................................................
    struct deriv_index_set {
      unsigned int tindex;
      unsigned int t_prime_index;
      unsigned int tau_index;
      INDEX npi;
      INDEX diff;
      INDEX ki;
      INDEX i;
      INDEX a;
      INDEX n;
      bool initialguess;
      bool regularized;
      deriv_index_set(unsigned int tindex, unsigned int t_prime_index, unsigned int tau_index, INDEX npi, 
          INDEX diff, INDEX ki, INDEX i, INDEX a, INDEX n, bool initialguess, bool regularized){
        this->tindex = tindex;
        this->t_prime_index = t_prime_index;
        this->tau_index = tau_index;
        this->npi = npi;
        this->diff = diff;
        this->ki = ki;
        this->i  = i;
        this->a  = a;
        this->n  = n;
        this->initialguess = initialguess;
        this->regularized    = regularized;
      }
    };
    typedef qgcomplex (ssystem::clocal::*sfunc)(const deriv_index_set &) const;
    // . ...................................................................................
    // data write/read functions 
    // . ...................................................................................
    // dump all data for the current tindex, and for all t' for unequal time quantities
    void dumpdata() 
    {
      if( dump_file_handles.L_out  != NULL ){
        fprintf( dump_file_handles.L_out, "=====================================\ndumping data for t=%u\n", tindex);
        this->fregdump( dump_file_handles.L_out, this->lop_register);
      } 
      for(unsigned int tprime = 0; tprime <= this->tindex; tprime++){
        if( dump_file_handles.R_out ){
          fprintf( dump_file_handles.R_out, "==========================================\n");
          fprintf( dump_file_handles.R_out, "dumping data for t=%u, t'=%u\n", tindex, tprime);
          fregdump( dump_file_handles.R_out, ssys.get_R_set( this->tindex, tprime ) );
        }
        if( dump_file_handles.S_out ){
          fprintf( dump_file_handles.S_out, "==========================================\n");
          fprintf( dump_file_handles.S_out, "dumping data for t=%u, t'=%u\n", tindex, tprime);
          fregdump( dump_file_handles.S_out, ssys.get_S_set( this->tindex, tprime ) );
        }
        if( dump_file_handles.PHI_out ){
          fprintf( dump_file_handles.PHI_out, "==========================================\n");
          fprintf( dump_file_handles.PHI_out, "dumping data for t=%u, t'=%u\n", tindex, tprime);
          fregdump( dump_file_handles.PHI_out, ssys.get_phi_set( this->tindex, tprime ) );
        }
        if( dump_file_handles.FPHI_out ){
          fprintf( dump_file_handles.FPHI_out, "==========================================\n");
          fprintf( dump_file_handles.FPHI_out, "dumping data for t=%u, t'=%u\n", tindex, tprime);
          fregdump( dump_file_handles.FPHI_out, ssys.get_F_phi_set( this->tindex, tprime ) );
        }
        if( dump_file_handles.G_out ){
          fprintf( dump_file_handles.G_out, "==========================================\n");
          fprintf( dump_file_handles.G_out, "dumping data for t=%u, t'=%u\n", tindex, tprime);
          fregdump( dump_file_handles.G_out, ssys.get_G_set( this->tindex, tprime ) );
        }
        if( dump_file_handles.FG_out ){
          fprintf( dump_file_handles.FG_out, "==========================================\n");
          fprintf( dump_file_handles.FG_out, "dumping data for t=%u, t'=%u\n", tindex, tprime);
          fregdump( dump_file_handles.FG_out, ssys.get_F_G_set( this->tindex, tprime ) );
        }
      }
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    inline qgcomplex get_RxPhi(const deriv_index_set & dset) const
    {
      const INDEX & npi = dset.npi;
      const INDEX & diff= dset.diff;
      const INDEX & ki  = dset.ki;
      const INDEX & i   = dset.i;
      const INDEX & a   = dset.a;
      const INDEX & n   = dset.n;
      const bool regularized = dset.regularized;
      const bool  & initialguess = dset.initialguess;
      //const bool  & equaltime    = dset.equaltime;
      unsigned int t_prime_index = dset.t_prime_index;
      unsigned int tindex        = dset.tindex;
      unsigned int tau_index1    = dset.tau_index; // used in R(t,t') as the 2nd argument
      unsigned int tau_index2    = dset.tau_index; // used in Phi(t,t') as the 1st argument
      const vector<int> & nkdom = (this->nkdomain);
      qgcomplex phival; 
      if( tau_index2 < t_prime_index){
        regset Phireg = this->phi_data->getvar( t_prime_index, tau_index2 ); 
        if(initialguess) Phireg.setprimary();
        else	         Phireg.setintermediate(); 
        int iindex = abs(static_cast<int>(diff));
        REAL factor = nkdom[iindex < static_cast<int>(nkdom.size()) ? iindex : 0]%2 == 0 ? 1 : -1;
        phival = factor*get_register(Phireg, diff, ki, n, a);
      }else{
        regset Phireg = this->phi_data->getvar( tau_index2, t_prime_index ); 
        if(initialguess) Phireg.setprimary();
        else	         Phireg.setintermediate(); 
        phival = get_register(Phireg, diff, ki, a, n); 
      }
      regset Rreg = regularized ? this->R_data_regularized->getvar(tindex, tau_index1) : 
        this->R_data->getvar( tindex, tau_index1);
      if(initialguess){
        Rreg.setprimary();
      }else{
        Rreg.setintermediate();
      }
      const qgcomplex & rval = get_register(Rreg, npi, ki, i, a);
      return rval*phival;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    inline qgcomplex get_RxG(const deriv_index_set & dset) const
    {
      // a few aliases:
      const INDEX & npi = dset.npi;
      const INDEX & diff= dset.diff;
      const INDEX & ki  = dset.ki;
      const INDEX & i   = dset.i;
      const INDEX & a   = dset.a;
      const INDEX & n   = dset.n;
      const bool  & initialguess = dset.initialguess;
      //const bool  & equaltime    = dset.equaltime;
      unsigned int t_prime_index = dset.t_prime_index;
      unsigned int tindex        = dset.tindex;
      unsigned int tau_index1    = dset.tau_index; // used in R(t,t') as the 2nd argument
      unsigned int tau_index2    = dset.tau_index; // used in G(t,t') as the 1st argument
      regset Greg = this->G_data->getvar(tau_index2, t_prime_index);
      regset Rreg = this->R_data->getvar( tindex, tau_index1);
      if(initialguess){
        Greg.setprimary();
        Rreg.setprimary();
      }else{
        Greg.setintermediate();
        Rreg.setintermediate();
      }
      const qgcomplex & Gval   = get_register(Greg, diff, ki, a, n);
      const qgcomplex & rval   = get_register(Rreg, npi, ki, i, a);
      return rval*Gval;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    inline qgcomplex get_SxG(const deriv_index_set & dset) const 
    { // a few aliases:
      const INDEX & npi = dset.npi;
      const INDEX & diff= dset.diff;
      const INDEX & ki  = dset.ki;
      const INDEX & i   = dset.i;
      const INDEX & a   = dset.a;
      const INDEX & n   = dset.n;
      const bool regularized = dset.regularized;
      const bool  & initialguess = dset.initialguess;
      //const bool  & equaltime    = dset.equaltime;
      unsigned int t_prime_index = dset.t_prime_index;
      unsigned int tindex        = dset.tindex;
      unsigned int tau_index1    = dset.tau_index; // used in S(t,t') as the 2nd argument
      unsigned int tau_index2    = dset.tau_index; // used in G(t,t') as the 2nd argument
      regset Greg = this->G_data->getvar( t_prime_index, tau_index2);
      regset Sreg = regularized ? 
        this->S_data_regularized->getvar( tindex, tau_index1) :
        this->S_data->getvar( tindex, tau_index1);
      if(initialguess){
        Greg.setprimary();
        Sreg.setprimary();
      }else{
        Greg.setintermediate();
        Sreg.setintermediate();
      }
      REAL factor = 1;
      if(static_cast<size_t>(abs(diff)) < nkdomain.size()){
        int idiff = nkdomain[abs(diff)];
        if(idiff % 2 != 0) factor = -1; 
      }
      //REAL factor = diff % 2 == 0 ? 1 : -1;
      const qgcomplex & sval   = get_register(Sreg, npi, ki, i, a);
      const qgcomplex & gval   = get_register(Greg, diff, ki, n , a);
      return factor*sval*gval;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    inline qgcomplex get_RxPhi_plus_SxG(const deriv_index_set & dset) const
    {
      return (get_RxPhi(dset) + get_SxG(dset));
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    inline qgcomplex & get_register_nonconst(regset & reg, const INDEX nk, const INDEX ki, const INDEX i, const INDEX j )
    {
      INDEX mmode = nk; // m_ang_mode;
      if(mmode < 0){
        mmode=(num_nks-1)-nk; // for negative modes
      }
      //return reg[ rank4( mmode, k_domain_index, i_field_index, n_field_index ) ]; // use rank indexing rather than direct linear access
      unsigned int index = ((static_cast<size_t>(mmode)*this->num_ks
            +static_cast<size_t>(ki))*this->numfields + static_cast<size_t>(i))*this->numfields 
        +static_cast<size_t>(j); // it is  faster to index the data set directly, rather than constructing a rank4 object*/
      return reg.data[ index ];
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    inline const qgcomplex & get_register(const regset & reg, const INDEX nk, const INDEX ki, const INDEX i, const INDEX j) 
      const
      {
        int mmode = nk; // m_ang_mode;
        if(mmode < 0){ // apply the reality constraint automatically
          mmode=(num_nks-1)-nk; // m_ang_mode; // for negative modes
        }
        if( static_cast<unsigned int>(abs(nk)) >= num_nks ){
          return null_complex;
        }
        //return reg[ rank4( mmode, k_domain_index, i_field_index, n_field_index ) ]; // use rank indexing rather than direct linear access
        unsigned int index = ((static_cast<size_t>(mmode)*this->num_ks
              +static_cast<size_t>(ki))*this->numfields + static_cast<size_t>(i))*this->numfields 
          +static_cast<size_t>(j); // it is  faster to index the data set directly, rather than constructing a rank4 object*/
        return reg.data[index];
      }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    inline void assign_reg(regset & reg, const qgcomplex & val, const INDEX ang_mode, const INDEX ki, 
        const INDEX i, const INDEX n  )
    {
      qgcomplex & element = this->get_register_nonconst(reg, ang_mode, ki, i ,n);
      if(ang_mode != 0){
        qgcomplex & neg_element = this->get_register_nonconst(reg, 0-ang_mode, ki, i, n);
        element = val;
        const int & nk = nkdomain[ static_cast<size_t>(abs(ang_mode)) ];
        REAL coeff = nk  % 2 == 0 ? 1.0 : -1.0;
        neg_element = qgcomplex(coeff*element.real(), -1.0*coeff*element.imag() );
      }else{
        element = val;
      }
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    const qgcomplex & get_register_noinline(const regset & reg, const INDEX m_ang_mode, const INDEX k_domain_index, const INDEX i_field_index, const INDEX j_field_index )  const 
    {
      int mmode = m_ang_mode;
      if(mmode < 0){
        mmode=((this->nkdomain).size()-1)-m_ang_mode; // for negative modes
      }
      if( static_cast<unsigned int>(abs(m_ang_mode)) >= num_nks ){
        return null_complex;
      }
      /* unsigned int index = ((static_cast<size_t>(mmode)*this->kdomain_centered->size()
         +static_cast<size_t>(k_domain_index))*this->numfields + static_cast<size_t>(i_field_index))*this->numfields
         +static_cast<size_t>(n_field_index); // it is  faster to index the data set directly, rather than constructing a rank4 object
         return reg[index];
         */  
      return reg[ rank4( mmode, k_domain_index, i_field_index, j_field_index ) ]; // use rank indexing rather than direct linear access
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    const static_container * getScoeffs() const{
      return &(mo.scoeffs);
    };
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    const static_container * getRcoeffs() const{
      return &(mo.rcoeffs);
    };
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  public:
    typedef const regset & (ssystem::clocal::*get_reg_set_func)(unsigned int, unsigned int) const;
    typedef const regset & (ssystem::clocal::*get_reg_set_single_time_func)(unsigned int) const;


    // . ...................................................................................
    // compilation/initialization functions
    // . ...................................................................................
    bool generateinitialconditions_phi_fourier_angular(initialdata & idata)
    {
      // iterate through all rule definitions and apply values
      for(size_t iter = 0; iter < idata.rules.size(); iter++){
        const rank4 & rpos = idata.rules[iter].rpos;
        bool i_index_free = rpos[2] < 0; 
        bool j_index_free = rpos[3] < 0;
        bool k_index_free = rpos[1] < 0;
        bool nkindex_free= rpos[0] < 0;
        for(INDEX npi = nkindex_free ? 0 : rpos[0]; npi <= (nkindex_free ? nkdomain.size()-1 : rpos[0]); ++npi){
          idata.index_mode = static_cast<int_type>(npi);
          idata.modeval = static_cast<int_type>(nkdomain[npi]);
          for(INDEX ki = k_index_free ? 0 : rpos[1]; ki <= (k_index_free ? num_ks - 1: rpos[0]); ki++){
            idata.index_k = static_cast<int_type>(ki+1);
            const REAL & Dk = dkdomain[ki];
            idata.kval = static_cast<float_type>(kdomain_centered[ki]);
            for(INDEX i = i_index_free ? 0 : rpos[2]; i <= (i_index_free ? numfields - 1 : rpos[2]); i++){
              idata.index_i = static_cast<int_type>(i+1);
              for(INDEX j = j_index_free ? 0 : rpos[2]; j <= (j_index_free ? numfields - 1 : rpos[2]); j++){
                idata.index_j = static_cast<int_type>(j+1);
                const cmplx_type & pval = idata.rules[iter].parser.Eval().GetComplex();
                qgcomplex phi_val( pval.real() * Dk, pval.imag() * Dk );
                this->assign_reg(phi_registers, phi_val, npi, ki, i, j);
              }}}}
      }
      return true;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    bool initializefields(spreprocdata & spdata)
    {
      initialdata & idata = *(spdata.getdata());
      // now iterate through the set of rules to generate the initial conditions
      // on Phi
      // first zero out the initial field in case no initial conditions were defined
      regset & phireg = this->phi_data->getvar(0,0);
      phireg.setprimary();
      phireg.zeromemory();
      phireg.setintermediate();
      phireg.zeromemory();

      try {
        // initialize the 0th tstep register for G
        this->initG(0, this->stepequaltimeonly);
        if( ! this->generateinitialconditions_phi_fourier_angular( idata) ){ 
          return false;
        } 
      }catch(const ParserError & e){
        fprintf(stderr, "%s\n", e.what());
        fprintf(stderr, "SSYSTEM INIT ERROR: error in expression evaluation.\n");
        return false;
      }
      // check consistency
      this->check_realizability = mo.get_global_param_by_name(gdb_str_IGNOREREALIZABILITYCONSTRAINTS).getiValue() == 0;
      if( check_realizability ){
        if(!this->check_field_realizability(0,true)){
          fprintf(stderr, "Model initialization failed because initial conditions violate symmetry properties of the covariance field.\n");
          return false;
        };
      }
      // now go through and copy initial data into the intermediate register of the 1st time step,
      // this is needed for the euler modified corrector step
      (*this->phi_registers.get_intermediate_data()) = (*this->phi_registers.get_primary_data());
      bool nonlinear = mo.get_global_param_by_name(gdb_str_NONLINEAR).getiValue() != 0;
      if(!nonlinear){
        mo.printfmsg("Warning: Model has nonlinear terms turned off.\n");
      }
      // calculate R and S
      if(nonlinear){
        this->calculateR(0,0, false);
        this->calculateR(0,0, true);
        this->calculateS(0,0, false);
        this->calculateS(0,0, true); 
      }
      return true;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  public:
    void setnumbertimesteps(const unsigned int & nsteps)
    {
#ifdef DEBUG
      assert(nsteps >= 0);
      this->timestepwidths.resize(nsteps, -1);
#else
      if(this->tindex > 0 && this->timestepwidths != NULL){
        REAL * temp = new REAL[this->tindex+1];
        for(size_t i = 0; i < this->tindex; i++) temp[i] = this->timestepwidths[i];
        this->timestepwidths = new REAL[nsteps+1];
        for(size_t i = 0; i < this->tindex; i++) this->timestepwidths[i] = temp[i];
        delete [] temp;
      }else{
        this->timestepwidths = new REAL[nsteps+1];
      }
      for(size_t i = this->tindex; i < nsteps+1; i++){
        this->timestepwidths[i] = -1;
      };
#endif
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    void settimestepwidth(const unsigned int & tindex, const REAL & tstep){
      assert(tindex <= this->tindex);
      assert(tstep > 0);
      this->timestepwidths[tindex] = tstep;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    void addtimeQL(const int tindex, const REAL & timestep, const bool initialguess)
    {
      // expand time variables to accomodate the next step
      if( this->stepequaltimeonly){
        if(this->timestepwidths[tindex] < 0){ // meaning the registers have not yet been allocated
          this->phi_data->addeqtime();
          this->F_phi_data->addeqtime();
        }
      }else if( this->tprimeintegrationlist ){
        const vector<int> & tlist = *(tprimeintegrationlist);
        phi_data->addtimes(tlist);
        G_data->addtimes(tlist);
        F_phi_data->addtimes( tlist );
        F_G_data->addtimes( tlist ); 
        this->initG(tindex+1, this->stepequaltimeonly); // only call this function after adding new time to G
      }else{
        this->addtime(tindex, timestep, initialguess);
        return;
      }
      this->settimestepwidth(tindex, timestep);
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    void addtime(const int tindex, const REAL & timestep, const bool initialguess)
    {
      // expand time variables to accomodate the next step 
      if(initialguess){
        if(this->timestepwidths[tindex] < 0){ 
          // meaning the registers have not yet been allocated
          // this function may be called with the future registers already allocated if 
          // the user requests the model to rewindtime via the rewindtime(nsteps) diagostic function
          this->F_phi_data->addtime();
          if( this->useregularization ) this->F_phi_data_regularized->addtime();
          this->F_G_data->addtime();
          if(this->stepequaltimeonly){
            this->phi_data->addeqtime();
            this->G_data  ->addeqtime();
          }else{
            phi_data->addtime();
            G_data->addtime();
            if(nonlinear){
              this->R_data->addtime();
              this->S_data->addtime();
              if(this->useregularization){
                this->R_data_regularized->addtime();
                this->S_data_regularized->addtime();
              }
            }
          }
        }
        this->initG(tindex+1, this->stepequaltimeonly); // only call this function after adding new time to G
        this->settimestepwidth(tindex, timestep);
      }
    } 
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    void create_linearop_map(spreprocdata & spdata, ParserX & P)
    {
      for(size_t op_iter = 0; op_iter < spdata.OPdata.size(); op_iter++){ // iterate over all the operators
        if(spdata.OPdata[op_iter].evaluate_op){ // if the operator is supposed to be decomposed then,
          for(size_t el_iter = 0; el_iter < spdata.OPdata[op_iter].expressions.size(); el_iter++){
            const rank3 & pos = spdata.OPdata[op_iter].expressions[el_iter].first;
            const token_list & opstr = spdata.OPdata[op_iter].expressions[el_iter].second;
            const INDEX & ang_mode = pos[0];
            const rank2 fpos( pos[1], pos[2] );
            assert( fpos[0] < static_cast<INDEX>(numfields) && fpos[0] >= 0 );
            assert( fpos[1] < static_cast<INDEX>(numfields) && fpos[1] >= 0 );
            if( ang_mode >= 0 ){
              assert(ang_mode < static_cast<INDEX>(nkdomain.size()) );
              this->linear_ops.lop_list.push_back(lopinfo::element());
              lopinfo::element & el = linear_ops.lop_list.back();
              el.assign( ang_mode, fpos, opstr, P );
            }else{
              for(INDEX mi = 0; mi < static_cast<INDEX>(num_nks); mi++){
                this->linear_ops.lop_list.push_back(lopinfo::element());
                lopinfo::element & el = linear_ops.lop_list.back();
                el.assign( mi, fpos, opstr, P );
              }
            }
          }
        }
      }
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    void create_dissipop_map(spreprocdata & spdata, ParserX & P)
    {
      const vector< pair<INDEX, token_list> > & sdata = spdata.dissip_data;
      for(size_t iter = 0; iter < sdata.size(); iter++){
        INDEX ipos = sdata[iter].first;
        assert( ipos < static_cast<INDEX>(numfields) && ipos >= 0 );
        this-> dissipation_ops.push_back( pair< INDEX, ParserX>( ipos, P ) ); 
        // hopefully the copy constructors work correctly ?? !! for the parser 
        ParserX & pset = this->dissipation_ops.back().second;
        const token_list & expr = sdata[iter].second;
        pset.SetExpr( expr );
        try{ 
          pset.Eval(); // force compilation ...
        }catch(const ParserError & e){
          fprintf(stderr, "ParserX::exception with message: %s.\n", e.what());
          fatal_stop("Parser evaluation failed at unexpected location.");
        }
      }
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    void create_inhomogeneousop_map(spreprocdata & spdata)
    {
      const vector< spreprocdata::ihomo_data > & sdata = spdata.inhomogeneous_data;
      for(size_t iter = 0; iter < sdata.size(); iter++){
        const rank3 & pos = sdata[iter].pos;
        const INDEX & mode_ref = pos[0];
        const INDEX & ipos_ref = pos[1];
        const INDEX & jpos_ref = pos[2];
        const string ang_mode_var_name = sdata[iter].ang_mode_var_name.size() > 0 ? sdata[iter].ang_mode_var_name : "m";
        const string k_var_name = sdata[iter].k_var_name.size() > 0 ? sdata[iter].k_var_name : "k";
        const string field_i_var_name = sdata[iter].field_i_var_name.size() > 0 ? sdata[iter].field_i_var_name : "i";
        const string field_j_var_name = sdata[iter].field_j_var_name.size() > 0 ? sdata[iter].field_j_var_name : "j";
        INDEX mode_min, mode_max, ipos_min, ipos_max, jpos_min, jpos_max;
        if(mode_ref < 0 ){
          mode_min = 0; mode_min = nkdomain.back(); 
        }else{ mode_min = mode_ref; mode_max = mode_ref; }
        if(ipos_ref < 0 ){
          ipos_min = 1; ipos_max = numfields;
        }else{ ipos_min = ipos_ref; ipos_max = ipos_ref; }
        if(jpos_ref < 0){
          jpos_min = 1; jpos_max = numfields;
        }else{ jpos_min = jpos_ref; jpos_max = jpos_ref; }
        for( INDEX mode = mode_min; mode <= mode_max; mode++){
          if( mode > this->nkdomain.back() ){ // out of range
            fprintf(stderr, "Angular mode '%i' is out of range.\n", static_cast<int>(mode));
            fprintf(stderr, "Error defining inhomogeneous operator '%s'.\n", sdata[iter].name.c_str());
            throw cgcmexception("S-system failed to initialize.", -1);
            return;
          }
          for( INDEX ipos = ipos_min; ipos <= ipos_max; ipos++){
            for( INDEX jpos = jpos_min; jpos <= jpos_max; jpos++){
              if( ipos > static_cast<INDEX>(numfields) || ipos < 1 ){
                fprintf(stderr, "Field index '%i' is out of range.\n", static_cast<int>(ipos));
                fprintf(stderr, "Error defining inhomogeneous operator '%s'.\n", sdata[iter].name.c_str());
                throw cgcmexception("S-system failed to initialize.", -1);
                return;
              }
              if( jpos > static_cast<INDEX>(numfields) || jpos < 1 ){
                fprintf(stderr, "Field index '%i' is out of range.\n", static_cast<int>(jpos));
                fprintf(stderr, "Error defining inhomogeneous operator '%s'.\n", sdata[iter].name.c_str());
                throw cgcmexception("S-system failed to initialize.", -1);
                return;
              }
              ParserX P(pckNONE);
              P.DefineVar(k_var_name, (this->reg.get_k()));
              P.DefineVar(ang_mode_var_name, (this->reg.get_m()));
              P.DefineVar(field_i_var_name, (this->reg.get_i()));
              P.DefineVar(field_j_var_name, (this->reg.get_j()));
              mo.assignparserdata(P, "c", true);
              mo.assigndiagfunctions(P);
              this->inhomogeneous_ops.push_back( pair<rank3, ParserX>( rank3(mode, ipos-1, jpos-1), P ) );
              ParserX & pset = this->inhomogeneous_ops.back().second;
              const token_list & expr = sdata[iter].expr;
              pset.SetExpr(expr);
              try{
                pset.Eval();
              }catch( const ParserError & e){
                fprintf(stderr, "Parser::exception with message: %s.\n", e.what());
                fatal_stop("Parser evaluation failed at unexpected location.");
              }
            }}
        }
      }
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    void initG(unsigned int tstep, bool equaltimeonly)
    {
      const unsigned int & nfs = numfields;
      const unsigned int & ksize = num_ks;
      const vector<REAL> & kdom = dkdomain;
      for(unsigned int tprime = ( equaltimeonly ? tstep : 0 ); tprime <= tstep; tprime++){
        this->G_registers = this->G_data->getvar(tstep, tprime);
        this->G_registers.setprimary();
        this->G_registers.zeromemory();
        this->G_registers.setintermediate();
        this->G_registers.zeromemory();
        if( tprime == tstep ){
          int nki = 0;
          for(unsigned int ki = 0; ki < ksize; ki++){
            // get the K width
            const REAL dk = kdom[ki];
            for(unsigned int i = 0; i < nfs; i++){
              this->G_registers.setprimary();
              this->assign_reg( this->G_registers, qgcomplex(dk,0), nki, ki, i, i);
              this->G_registers.setintermediate();
              this->assign_reg( this->G_registers, qgcomplex(dk,0), nki, ki, i, i);
            }  
          }
        }
      }
    }
    //......................................................................................
    // integration functions:
    //......................................................................................
    qgcomplex trap_rule_integrate(sfunc function, deriv_index_set & dset, unsigned int t0, unsigned int tmax)
    {
      qgcomplex rvalue(0.0,0.0);
      for(unsigned int tau_index = t0; tau_index < tmax; tau_index++){
        dset.tau_index = tau_index;
        qgcomplex fn = (*this.*function)(dset);
        dset.tau_index = tau_index+1;
        qgcomplex fn_next = (*this.*function)(dset);
        const REAL timestep = (this->timestepwidths[tau_index]*0.5);
        assert(timestep > 0);
        rvalue += qgcomplex(timestep,0)*(fn+fn_next);
        check_for_signals(); // check for INTERRUPT signals
      }
      return rvalue;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    bool check_field_realizability(unsigned int t1, bool displaymessage) const{
      // also prints warnings if the mirror symmetry is not present in the initial conditions
      // check that Phi_{ij}^{m} = (-1)^m Phi_{ji}^{m} (Property of Phi) for the equal time only (unequal does not exist for t'>t
      unsigned int t2 = t1; // equal time only
      regset Phireg1 = this->phi_data->getvar(t1, t1); 
      regset Phireg2 = this->phi_data->getvar(t1, t1);
      const vector<int> & nkdom = (this->nkdomain);
      const vector<REAL> & kdom = (this->kdomain_centered);
      Phireg1.setprimary();
      Phireg2.setprimary(); 
      for(int loop = 0; loop < 2; loop++){ // check both primary and intermediate (provisional values)
        if(loop == 1){ 
          Phireg1.setintermediate();
          Phireg2.setintermediate();
        }
        for(INDEX ki = 0; ki < static_cast<INDEX>(this->num_ks); ki++){
          for(INDEX nki = 0; nki < static_cast<INDEX>(num_nks); nki++){ 
            for(INDEX i = 0; i < static_cast<INDEX>(numfields); i++){
              for(INDEX j = 0; j < static_cast<INDEX>(numfields); j++){
                const qgcomplex & val1 = this->get_register(Phireg1, nki, ki, i, j);
                const qgcomplex & val2 = this->get_register(Phireg2, nki, ki, j, i);
                if(val1 != val1) continue; // nan values, ignore
                if(val2 != val2) continue; // nan values ignore
                REAL coeff = nkdom[nki] % 2 == 0 ? 1.0 : -1.0;
                bool normalizer_zero = abs(val1.real()) < 0.000000001;
                REAL rdiff = ! normalizer_zero ?
                  abs((val1.real() - coeff*val2.real())/val1.real()) :
                  abs((val1.real() - coeff*val2.real()));
                REAL idiff = ! normalizer_zero ?
                  abs((val1.imag() - coeff*val2.imag())/val1.imag()) :
                  abs((val1.imag() - coeff*val2.imag()));
                if( rdiff > 0.00000001 || idiff > 0.00000001 ){ 
                  if(displaymessage){
                    fprintf(stderr, "Error: Relative diff real = %f, relative diff imag = %f.\n", rdiff, idiff );
                    fprintf(stderr, "Error: Phi[nk=%i, k=%f, i=%i, j=%i] (t=%u, t'=%u) := (%f, %f) ",
                        nkdom[nki], kdom[ki], i+1, j+1, t1, t2, val1.real(), val1.imag());
                    fprintf(stderr," !=  (-1)^(nk) Phi[nk=%i, k=%f, i=%i, j=%i] (t'=%u, t=%u) := (%f, %f) \n",
                        nkdom[nki], kdom[ki], j+1, i+1, t2, t1, val2.real(), val2.imag());
                  }
                  return false;
                }
                if(displaymessage){
                  int nk = nkdom[nki];
                  if(nk % 2 == 0){ // make sure there is no imag part
                    if(val1.imag() != 0.0){
                      fprintf(stderr, "Warning: Phi[nk=%i, k=%f, i=%i, j=%i] (t=%u, t'=%u) is not obeying the mirror symmetry.\n",
                          nk, kdom[ki], i+1, j+1, t1, t2);
                    }
                  }else{ // make sure there is no real part
                    if(val1.real() != 0.0){
                      fprintf(stderr, "Warning: Phi[nk=%i, k=%f, i=%i, j=%i] (t=%u, t'=%u) is not obeying the mirror symmetry.\n",
                          nk, kdom[ki], i+1, j+1, t1, t2);
                    }
                  }
                }  
              }}
          }
        }
      }
      return true;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    REAL getflux(unsigned int index_i, unsigned int index_j, unsigned int tindex, bool initial_guess) const
    {
      const REAL coeff = fourier_transform_coefficient*(fourier_transform_coefficient);
      const vector<REAL> & k2 = this->kdomain_centered_squared; 
      REAL area = 0; 
      int kc_index = 0;
      // the phis that are evolved in the model are actually Phi(k)*dk
      INDEX mode = -1;
      const vector<int> & nkdom = (this->nkdomain);
      for(size_t iter = 0; iter < nkdom.size(); iter++){
        if( nkdom[iter] == 1 ){
          mode = static_cast<INDEX>(iter);
          break; 
        }
      }
      if( mode < 0 ){
        return 0.0;
      };
      regset phi;
      if(initial_guess){
        phi = this->phi_data->getvar(tindex+1, tindex+1);
        phi.setintermediate();
      }else{
        phi = this->phi_data->getvar(tindex, tindex);
        phi.setprimary();
      }
      for(unsigned int kindex = 0; kindex < k2.size(); kindex++){
        const qgcomplex & fval = get_register(phi, mode, kc_index, index_i, index_j);
        area+= k2[kindex]*fval.imag();
        kc_index++;
      }
      area *= coeff;
      return static_cast<double>(area);
    }
    bool init(spreprocdata & spdata)
    {
      /* ..................................................................................... */
      // we need space for both the positive and negative angular fourier modes
      // that's why the size of the 1st dimension is 2*nkdomain.size()
      rank4 edges( 2*nkdomain.size(), kdomain_centered.size(), numfields, numfields );

      this->phi_registers.allocate(edges); // allocate the first registers for t,t'=0, subsequent ones are allocated as needed
      // we call allocatefull to allocate both the primary and intermediate register set
      // subsequent allocations allocate only an additional primary register
      this->G_registers.allocate(edges);


      if(nonlinear){
        this->S_registers.allocate(edges); 
        this->R_registers.allocate(edges);
        if(this->useregularization){
          this->S_registers_regularized.allocate(edges);
          this->R_registers_regularized.allocate(edges);
        }
      }
      this->F_phi_registers.allocate(edges);
      if(this->useregularization) this->F_phi_registers_regularized.allocate(edges);

      this->F_G_registers.allocate(edges);

      this->lop_register.allocate(edges);
      this->lop_non_inverted.allocate(edges);

      this->dissip_register.allocate(edges);
      this->inhomo_forcing_register.allocate(edges);
      /* ..................................................................................... */
      // Assigning namelist parser data to parser template as well as the special variables k2,kx,ky
      this->reg.set_reg(kdomain_centered[0]); // initialize registers
      this->reg.set_ang(nkdomain[0]);
      this->reg.set_fields(1,1);


      ParserX P(pckNONE);
      P.DefineVar("k", (this->reg.get_k()));
      P.DefineVar("m", (this->reg.get_m()));
      P.DefineVar("i", (this->reg.get_i()));
      P.DefineVar("j", (this->reg.get_j()));
      mo.assignparserdata(P,"c",true);
      mo.assigndiagfunctions(P);
      this->create_linearop_map(spdata, P);

      ParserX Pv;
      Pv.DefineVar("k", (this->reg.get_k()));
      Pv.DefineVar("i", (this->reg.get_i()));
      mo.assignparserdata(Pv,"f",true);
      mo.assigndiagfunctions(Pv);
      this->create_dissipop_map(spdata, Pv);

      this->create_inhomogeneousop_map(spdata);

      unsigned int numthreads = 1;
#ifndef SERIALQGCM
      int NUMTHREADS = mo.get_num_threads();
      omp_set_num_threads(NUMTHREADS);
      numthreads = NUMTHREADS;
#endif

      // create two timed variable containers to contain the multi-time register sets 
      this->phi_data = new twotimevar(this->phi_registers);
      this->G_data = new twotimevar(this->G_registers);
      this->F_phi_data = new twotimevar(this->F_phi_registers);
      this->F_G_data = new twotimevar(this->F_G_registers);
      if(this->useregularization){
        this->F_phi_data_regularized = new twotimevar(this->F_phi_registers_regularized);
      }
      if(nonlinear){
        this->S_data = new twotimevar(this->S_registers);
        this->R_data = new twotimevar(this->R_registers);
        if(this->useregularization){
          this->R_data_regularized = new twotimevar(this->R_registers_regularized);
          this->S_data_regularized = new twotimevar(this->S_registers_regularized);
        }
      }

      this->reg.set_reg((kdomain_centered)[0]); // init the system registers
      if(!this->initializefields( spdata)) return false;

      vector< pair<const char *, FILE ** > > dplist;
      dplist.push_back( pair<const char *, FILE**>( gdb_str_L_DUMP, &(this->dump_file_handles.L_out) ) );
      dplist.push_back( pair<const char *, FILE**>( gdb_str_FPHI_DUMP, &(this->dump_file_handles.FPHI_out) ) );
      dplist.push_back( pair<const char *, FILE**>( gdb_str_PHI_DUMP, &(this->dump_file_handles.PHI_out) ) );
      dplist.push_back( pair<const char *, FILE**>( gdb_str_R_DUMP, &(this->dump_file_handles.R_out) ) );
      dplist.push_back( pair<const char *, FILE**>( gdb_str_S_DUMP, &(this->dump_file_handles.S_out) ) );
      dplist.push_back( pair<const char *, FILE**>( gdb_str_G_DUMP, &(this->dump_file_handles.G_out) ) );
      dplist.push_back( pair<const char *, FILE**>( gdb_str_FG_DUMP, &(this->dump_file_handles.FG_out) ) );
      for(size_t it = 0; it < dplist.size(); it++){
        *(dplist[it].second) = NULL;
        if( mo.does_param_exist( dplist[it].first ) ){
          const parameter & param = mo.get_global_param_by_name(dplist[it].first );
          if(param.getType() == parameter::String){
            FILE * new_file = fopen( mo.get_global_param_by_name( dplist[it].first ).getsValue(), "w" );
            if( new_file ){
              *(dplist[it].second) = new_file;
            }else{
              fprintf(stderr, "Could not open file '%s' specified as an output file.\n", dplist[it].first );
              throw cgcmexception("S-system failed to initialize.", -1);
              return false;
            }
          }
        }
      }





      this->initialized = true;
      return true;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    void compute_lops(bool initialguess)
    {
      if(initialguess){
        this->lop_register.setprimary();
        this->lop_non_inverted.setprimary();
      }else{
        this->lop_register.setintermediate();
        this->lop_non_inverted.setintermediate();
      } 
      this->lop_register.zeromemory();
      this->lop_non_inverted.zeromemory(); 
      if( linear_ops.lop_list.size() <= 0 ) return;
      // iterate through the LOP info list
      // the calculation is a speckle sum (add to those elements that are nonzero)
      for(size_t iter = 0; iter < linear_ops.lop_list.size(); iter++){
        const INDEX & m = linear_ops.lop_list[iter].angular_mode;
        const INDEX & i = linear_ops.lop_list[iter].field_indices[0];
        const INDEX & j = linear_ops.lop_list[iter].field_indices[1];
        this->reg.set_ang(nkdomain[m]);
        this->reg.set_fields(i, j);
        for(INDEX ki = 0; ki < static_cast<INDEX>(num_ks); ki++){
          this->reg.set_reg(kdomain_centered[ki]);
          double rval = linear_ops.lop_list[iter].real();
          double ival = linear_ops.lop_list[iter].imag();
          qgcomplex new_val( static_cast<REAL>(rval), static_cast<REAL>(ival));
          new_val += this->get_register(this->lop_non_inverted, m, ki, i, j);
          this->assign_reg(this->lop_non_inverted, new_val, m, ki, i, j);
        };
      };
      // now we need to invert the lop using the Qinverses, which is isotropic
      const map_container<REAL, 3> & QI = mo.QIcoeffs;
      for(INDEX mi = 0; mi < static_cast<INDEX>( num_nks ); mi ++ ){
        for(INDEX ki = 0; ki < static_cast<INDEX>( this->num_ks ); ki++){
          for(INDEX i = 0; i < static_cast<INDEX>(numfields); i++){ 
            for(INDEX n = 0; n < static_cast<INDEX>(numfields); n++){
              qgcomplex sum(0,0);
              for(INDEX a = 0; a < static_cast<INDEX>(numfields); a++){
                const REAL & QIval = QI.getvalue(rank3(ki, i, a));
                const qgcomplex & lop_val = this->get_register(this->lop_non_inverted, mi, ki, a, n);
                sum+= QIval*lop_val; 
              }
              this->assign_reg(this->lop_register, sum, mi, ki, i, n); 
            }
          }
        } 
      }
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    void compute_dissip(bool initialguess)
    {
      if(initialguess){
        this->dissip_register.setprimary();
      }else{
        this->dissip_register.setintermediate();
      }
      this->dissip_register.zeromemory();
      for(INDEX ki = 0; ki < static_cast<INDEX>(num_ks); ki++){
        for(size_t iter = 0; iter < dissipation_ops.size(); iter++){
          const INDEX & i = dissipation_ops[iter].first; 
          ParserX & p = dissipation_ops[iter].second;
          this->reg.set_fields(i, i);
          this->reg.set_reg(kdomain_centered[ki]);
          const qgcomplex & val = p.Eval().GetComplex();
          qgcomplex new_val = this->get_register(this->dissip_register, 0, ki, i, i) + val;
          this->assign_reg(this->dissip_register, new_val, 0, ki, i, i); 
        }
      }
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    void compute_inhomogeneous_part(bool initial_guess)
    {
      if(initial_guess){
        this->inhomo_forcing_register.setprimary();
      }else{
        this->inhomo_forcing_register.setintermediate();
      }
      this->inhomo_forcing_register.zeromemory(); // speckle filling
      for(INDEX ki = 0; ki < static_cast<INDEX>(num_ks); ki++){
        for(size_t iter = 0; iter < inhomogeneous_ops.size(); iter++){
          const rank3 & pos = inhomogeneous_ops[iter].first; 
          ParserX & p = inhomogeneous_ops[iter].second;
          INDEX mode = pos[0];
          INDEX i = pos[1];
          INDEX j = pos[2];
          this->reg.set_ang( mode );
          this->reg.set_fields(i, j);
          this->reg.set_reg(kdomain_centered[ki]);
          const qgcomplex & val = p.Eval().GetComplex();
          qgcomplex new_val = this->get_register(this->inhomo_forcing_register, mode, ki, i, j) + val;
          this->assign_reg(this->inhomo_forcing_register, new_val, mode, ki, i, j); 
        }
      }
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    void stepforwardNL(REAL timestep, bool initialguess)
    {
      this->phi_registers = this->phi_data->getvar(tindex, tindex); // set reg to current time single time phi
      this->G_registers   = this->G_data  ->getvar(tindex, tindex);
      // expand time variables to accomodate the next step
      this->addtime(tindex, timestep, initialguess);
      // expands the time if we are doing a provisional step (initialguess)
      dbgprintf("starting new time integration %i, initialguess='%s'\n", this->tindex, initialguess ? "true" : "false");
      // check realizability
      {
        unsigned int t1 = tindex;
        //  if( check_realizability && ! this->check_field_realizability(t1, false) ){
        if( check_realizability && ! this->check_field_realizability(t1, true) ){
          fatal_stop("Phi field failed symmetry test during integration.");
        }
      }
      unsigned int cindex = initialguess ? tindex : tindex+1; // provisional stuff is stored in the intermediate register of the
      // next time step 
      unsigned int next_t_index = tindex + 1;
      // compute the angular mode Lij components for the current time step (depend only on tindex)
      this->compute_lops(initialguess);
      this->compute_dissip(initialguess);
      this->compute_inhomogeneous_part(initialguess); 
      // calculate the unequal time derivatives ( F(t,t') for all t' < t ) 
      {
#ifndef SERIALQGCM
        dbgprintf("starting parallel section.\n"); 
#pragma omp parallel 
#pragma omp for schedule(dynamic, 1)
#endif
        for(unsigned int t_prime_index = 0; t_prime_index <= tindex; t_prime_index++){
          try{ 
            regset Fphi = this->F_phi_data->getvar(cindex, t_prime_index);
            regset FG   = this->F_G_data->getvar(cindex, t_prime_index);
            //int tindex_dummy = initialguess ? tindex : tindex + 1; 
            this->calculate_F_phi(cindex, t_prime_index, initialguess, nonlinear, false, Fphi);
            this->calculate_F_G  (cindex, t_prime_index, initialguess, nonlinear, FG);
            // the last_Phi is always the corrected (primary) value at tindex (not cindex)
            regset last_Phi = this->phi_data->getvar(tindex, t_prime_index);
            regset last_G   = this->G_data->getvar(tindex, t_prime_index);
            regset next_Phi = this->phi_data->getvar( next_t_index, t_prime_index); 
            regset next_G   = this->G_data->getvar( next_t_index, t_prime_index);
            if(initialguess){
              advance_and_assign(initialguess, last_Phi, next_Phi, Fphi, Fphi);
              advance_and_assign(initialguess, last_G , next_G, FG, FG);
            }else{
              regset Fphi2 = this->F_phi_data->getvar(tindex, t_prime_index);
              regset FG2   = this->F_G_data->getvar(tindex, t_prime_index);
              advance_and_assign(initialguess, last_Phi, next_Phi, Fphi2, Fphi);
              advance_and_assign(initialguess, last_G , next_G, FG2, FG);
            }
#ifndef SERIALQGCM
            dbgprintf("thread '%i' ", omp_get_thread_num());
#endif
            dbgprintf("calculating R & S(%i, %i) for tindex = %i, initialguess = '%s'\n", 
                next_t_index, t_prime_index, this->tindex,
                initialguess ? "true" : "false");
            if(nonlinear){
              this->calculateR(next_t_index, t_prime_index, initialguess);
              this->calculateS(next_t_index, t_prime_index, initialguess);
            };
          }catch(const cgcmexception & ex){
            if( is_exception_SIG_INT(ex) ){ 
#ifdef SERIALQGCM
              break ;
#endif 
            }else{ throw ex; };
          }
        }
#ifndef SERIALQGCM
        dbgprintf("Ending parallel section.\n"); 
#endif
      } 
      // this section is now done serially

      check_for_signals(); // check for interrupt signals

      // now we do the equal time, G equal times do not evolve
      {
        regset Fphi = this->F_phi_data->getvar(cindex, cindex);
        regset Fphi_regularized;
        if(this->useregularization){
          Fphi_regularized = this->F_phi_data_regularized->getvar(cindex, cindex);
          this->calculate_F_phi(cindex, cindex, initialguess, nonlinear, true, Fphi_regularized);
        }
        // regset FG   = this->F_G_data->getvar(cindex, cindex);
        //int tindex_dummy = initialguess ? tindex : tindex + 1; 
        if(!initialguess){ // we have not yet calculated these, so do it
          this->calculate_F_phi(cindex, cindex, initialguess, nonlinear, false, Fphi);
        }
        // the last_Phi is always the corrected (primary) value at tindex (not cindex)
        regset last_Phi = this->phi_data->getvar(tindex, tindex);

        regset next_Phi = this->phi_data->getvar( next_t_index, next_t_index); 
        if(initialguess){
          advance_and_assign_equal_time( initialguess, last_Phi, next_Phi, Fphi_regularized, Fphi_regularized, Fphi, Fphi );
        }else{
          regset Fphi2 = this->F_phi_data->getvar(tindex, tindex);
          regset Fphi2_regularized;
          if(this->useregularization) Fphi2_regularized = this->F_phi_data_regularized->getvar(tindex, tindex);
          advance_and_assign_equal_time( initialguess, last_Phi, next_Phi, Fphi2_regularized, Fphi_regularized, Fphi2, Fphi);
        }
        if(nonlinear){ // calculate all the Rs and Ss for the current tindex
          this->calculateR(next_t_index, next_t_index, initialguess);
          this->calculateS(next_t_index, next_t_index, initialguess);
        }
      }

#ifdef DEBUG
      if( !initialguess ) this->dumpdata();
#endif

      // make sure the register order is correct here
      // the s fluxes will be calculated next, if this is an initial_guess , the flux should be calculated
      // using values from the current intermediate register set
      this->phi_registers = this->phi_data->getvar(next_t_index, next_t_index);
      this->G_registers = this->G_data->getvar(next_t_index, next_t_index);
      if(!initialguess){ 
        this->tindex++;
        this->phi_registers.setprimary();
        this->G_registers.setprimary();
      }else{
        this->phi_registers.setintermediate();
        this->G_registers.setintermediate();
      }
      }
      // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      void stepforwardQL(REAL timestep, bool initialguess)
      {
        this->phi_registers = this->phi_data->getvar(tindex, tindex); // set reg to current time single time phi
        if( !this->stepequaltimeonly) this->G_registers   = this->G_data  ->getvar(tindex, tindex);
        // expand time variables to accomodate the next step
        this->addtimeQL(tindex, timestep, initialguess);
        // expands the time if we are doing a provisional step (initialguess)
        dbgprintf("starting new time integration %i, initialguess='%s'\n", this->tindex, initialguess ? "true" : "false");
        // check realizability
        {
          unsigned int t1 = tindex;
          //  if( check_realizability && ! this->check_field_realizability(t1, false) ){
          if( check_realizability && ! this->check_field_realizability(t1, true) ){
            fatal_stop("Phi field failed symmetry test during integration.");
          }
        }
        unsigned int cindex = initialguess ? tindex : tindex+1; // provisional stuff is stored in the intermediate register of the
        // next time step 
        unsigned int next_t_index = tindex + 1;
        // compute the angular mode Lij components for the current time step (depend only on tindex)
        this->compute_lops(initialguess);
        this->compute_dissip(initialguess);
        this->compute_inhomogeneous_part(initialguess); 
        // calculate the unequal time derivatives ( F(t,t') for all t' < t ) 
        if(!this->stepequaltimeonly){
          // use integration list if available
          // for(unsigned int t_prime_index = 0; t_prime_index <= tindex; t_prime_index++){
          unsigned int t_prime_index = 0;
          unsigned int list_iter = 0;
          while( true ) {
            // first check loop condition (incrementation is done at the end to mimic a for loop
            if( this->tprimeintegrationlist ){
              if( list_iter >= tprimeintegrationlist->size() ) break;
              t_prime_index = (* tprimeintegrationlist)[list_iter];
            }
            if( t_prime_index > tindex ) break;
            try{ 
              regset Fphi = this->F_phi_data->getvar(cindex, t_prime_index);
              regset FG   = this->F_G_data->getvar(cindex, t_prime_index);
              //int tindex_dummy = initialguess ? tindex : tindex + 1; 
              this->calculate_F_phi(cindex, t_prime_index, initialguess, nonlinear, false, Fphi);
              this->calculate_F_G  (cindex, t_prime_index, initialguess, nonlinear, FG);
              // the last_Phi is always the corrected (primary) value at tindex (not cindex)
              regset last_Phi = this->phi_data->getvar(tindex, t_prime_index);
              regset last_G   = this->G_data->getvar(tindex, t_prime_index);
              regset next_Phi = this->phi_data->getvar( next_t_index, t_prime_index); 
              regset next_G   = this->G_data->getvar( next_t_index, t_prime_index);
              if(initialguess){
                advance_and_assign(initialguess, last_Phi, next_Phi, Fphi, Fphi);
                advance_and_assign(initialguess, last_G , next_G, FG, FG);
              }else{
                regset Fphi2 = this->F_phi_data->getvar(tindex, t_prime_index);
                regset FG2   = this->F_G_data->getvar(tindex, t_prime_index);
                advance_and_assign(initialguess, last_Phi, next_Phi, Fphi2, Fphi);
                advance_and_assign(initialguess, last_G , next_G, FG2, FG);
              }
              dbgprintf("calculating R & S(%i, %i) for tindex = %i, initialguess = '%s'\n", 
                  next_t_index, t_prime_index, this->tindex,
                  initialguess ? "true" : "false");
            }catch(const cgcmexception & ex){
              if( is_exception_SIG_INT(ex) ){ 
                break ;
              }else{ throw ex; };
            }
            // select next t' index
            if( this->tprimeintegrationlist ){
              list_iter ++;
            }else{
              t_prime_index ++;
            } 
            // end selection
          }
        } // end of if(stepequaltimeonly) 
        // this section is now done serially

        check_for_signals(); // check for interrupt signals

        // now we do the equal time, G equal times do not evolve
        {
          regset Fphi = this->F_phi_data->getvar(cindex, cindex);
          regset Fphi_regularized;
          if(this->useregularization){
            Fphi_regularized = this->F_phi_data_regularized->getvar(cindex, cindex);
            this->calculate_F_phi(cindex, cindex, initialguess, nonlinear, true, Fphi_regularized);
          }
          // regset FG   = this->F_G_data->getvar(cindex, cindex);
          //int tindex_dummy = initialguess ? tindex : tindex + 1; 
          if(!initialguess){ // we have not yet calculated these, so do it
            this->calculate_F_phi(cindex, cindex, initialguess, nonlinear, false, Fphi);
          }
          // the last_Phi is always the corrected (primary) value at tindex (not cindex)
          regset last_Phi = this->phi_data->getvar(tindex, tindex);

          regset next_Phi = this->phi_data->getvar( next_t_index, next_t_index); 
          if(initialguess){
            advance_and_assign_equal_time( initialguess, last_Phi, next_Phi, Fphi_regularized, Fphi_regularized, Fphi, Fphi );
          }else{
            regset Fphi2 = this->F_phi_data->getvar(tindex, tindex);
            regset Fphi2_regularized;
            if(this->useregularization) Fphi2_regularized = this->F_phi_data_regularized->getvar(tindex, tindex);
            advance_and_assign_equal_time( initialguess, last_Phi, next_Phi, Fphi2_regularized, Fphi_regularized, Fphi2, Fphi);
          }
        }

#ifdef DEBUG
        if( !initialguess ) this->dumpdata();
#endif

        // make sure the register order is correct here
        // the s fluxes will be calculated next, if this is an initial_guess , the flux should be calculated
        // using values from the current intermediate register set

        this->phi_registers = this->phi_data->getvar(next_t_index, next_t_index);
        if( ! this->stepequaltimeonly ){
          this->G_registers = this->G_data->getvar(next_t_index, next_t_index);
        }
        if(!initialguess){ 
          this->tindex++;
          this->phi_registers.setprimary();
          if(! this->stepequaltimeonly) this->G_registers.setprimary();
        }else{
          this->phi_registers.setintermediate();
          if(! this->stepequaltimeonly) this->G_registers.setintermediate();
        }
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void stepforward(REAL timestep, bool initialguess)
        {
          if(nonlinear){
            this->stepforwardNL(timestep, initialguess);
          }else{
            this->stepforwardQL(timestep, initialguess);
          }
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void int_dynamic_pq(
            const static_container & static_coeffs, regset & foutput, 
            regset & foutput_regularized, 
            bool R, const regset & Pregset, const regset & Qregset) 
        {
          // this function may be called in a parallel shared memory environment
          // in which case everything except writes to the destination regset (output) must be atomic
          // care must be taken that the destination is not read from or written to by other threads
          const rank11 * index_map =  static_coeffs.get_index_map();
          const int size = static_coeffs.get_size();
          foutput.zeromemory(); // set all values to 0 because this is a speckle integration
          foutput_regularized.zeromemory();
          regset output = foutput;
          const qgcomplex * static_coeff_data = static_coeffs.getcoefficients();
          for(int iter = 0; iter < size; iter++){
            //check_for_signals(); // check for interrupts
            const rank11 & r11 = index_map[iter];
            const INDEX & mki = r11[0]; //.mk;
            const INDEX & npi = r11[1]; //.nk;
            const INDEX & ki = r11[2]; //.ki;
            const INDEX & pi = r11[3]; //.pi;
            const INDEX & qi = r11[4]; //.qi;
            const INDEX &  i = r11[5];
            const INDEX &  j = r11[6]; 
            const INDEX &  k = r11[7]; 
            const INDEX &  m = r11[8]; 
            const INDEX &  a = r11[9]; 
            const INDEX &  b = r11[10];
            // check regularization condition
            if(this->useregularization){
              const REAL kmin_index = this->kminimum[static_cast<int>(ki)];
              if(pi < kmin_index && qi < kmin_index){
                output = foutput_regularized;
              }else{
                output = foutput;
              }
            }
            // You cannot write to pieces of a register directly because you will not write the negative modes!
            // Please use function assign_reg to assign the value to a register correctly ...
            const qgcomplex & element = R ? this->get_register( output, mki, ki, i, a ) :
              this->get_register( output, mki, ki, i, m );
            const qgcomplex & static_coeff = static_coeff_data[iter];
            INDEX diff =  mki - npi ;
            const qgcomplex & Qreg = this->get_register(Qregset, diff, qi, k, b ); 
            const qgcomplex & Preg = R ? this->get_register( Pregset, npi, pi, j, m ) :
              this->get_register( Pregset, npi, pi, j, a);
            const qgcomplex nval = static_coeff*Qreg*Preg;
            if(R){
              this->assign_reg( output, element + nval, mki, ki, i, a );
            }else{
              this->assign_reg( output, element + nval, mki, ki, i, m );
            }
          }
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void calculateR(unsigned int tindex, unsigned int t_prime_index, bool provisional)
        {
          // these functions must be atomic (thread safe) !!! with the following qualification:
          // the data written to Sreg must NOT be accessed by other threads
          // when using this function in parallel mode
          regset Rreg_regularized;
          if(this->useregularization) Rreg_regularized = this->R_data_regularized->getvar(tindex, t_prime_index);
          regset Rreg = this->R_data->getvar(tindex, t_prime_index);
          regset Greg = this->G_data->getvar(tindex, t_prime_index);
          regset Phireg=this->phi_data->getvar(tindex, t_prime_index);
          if(!provisional){  
            Rreg.setprimary();  Rreg_regularized.setprimary();
            Greg.setprimary();
            Phireg.setprimary();
          }else{
            Rreg.setintermediate(); Rreg_regularized.setintermediate();
            Greg.setintermediate();
            Phireg.setintermediate();
          }
          //dbgprintf("\n---- Beginning R calculation ----\n");
          int_dynamic_pq(*(this->getRcoeffs()), Rreg, Rreg_regularized, true, Greg, Phireg );
          //dbgprintf("---- Finished R calculation ----\n\n");
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void calculateS(unsigned int tindex, unsigned int t_prime_index, bool provisional)
        {
          // these functions must be atomic (thread safe) !!! with the following qualification
          // note that the data written to Sreg must NOT be accessed by other threads
          // when using this function in parallel mode
          regset Sreg_regularized;
          if(this->useregularization) Sreg_regularized = this->S_data_regularized->getvar(tindex, t_prime_index);
          regset Sreg = this->S_data->getvar(tindex, t_prime_index);
          regset Phireg=this->phi_data->getvar(tindex, t_prime_index);
          if(!provisional){
            Sreg.setprimary();  Sreg_regularized.setprimary();
            Phireg.setprimary();
          }else{
            Sreg.setintermediate(); Sreg_regularized.setintermediate();
            Phireg.setintermediate();
          }
          //dbgprintf("\n---- Beginning S calculation ----\n");
          int_dynamic_pq(*(this-> getScoeffs()), Sreg, Sreg_regularized, false, Phireg, Phireg );
          //dbgprintf("---- Finished R calculation ----\n\n");
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void calculate_F_phi(unsigned int tindex, unsigned int t_prime_index, bool initialguess, bool nonlinear, 
            bool regularized, regset Fdest){
          if(initialguess){ // assign to primary if the derivative is to be calculated is the for the provisional timestep
            // i.e. we are calculating the derivative at the corrected location ....
            Fdest.setprimary();
          }else{
            Fdest.setintermediate();
          }
          // Rs, Ss needed should already be calculated by this point
          INDEX nfs = static_cast<INDEX>(this->numfields);
          int nksize = num_nks; //static_cast<int>( this->nkdomain->size() );
          unsigned int ksize  = num_ks; //static_cast<unsigned int>( this->kdomain_centered->size() );
          // outer indices
          for(int nki = 0; nki < nksize; nki++){
            for(unsigned int ki = 0; ki < ksize; ki++){
              for(INDEX i = 0; i < nfs; i++){
                for(INDEX n = 0; n < nfs; n++){
                  // - - - - - iner contractions (summed over) index loops
                  qgcomplex sum(0,0);
                  for(INDEX a = 0; a < nfs; a++){
                    // do the linear part first ...
                    if(!regularized){ // we do not need to do the linear calculation twice!
                      this->calculate_F_linear_part(tindex, t_prime_index, sum, nki, ki, i, n, a, this->phi_data, initialguess);
                    }
                    // linear operations scope
                    // now do nonlinear part
                    if(nonlinear){
                      for(int mki = -nksize + 1; mki < nksize; mki++){
                        // we split into two pieces for simplicity tau = 0 to t' and tau = t' to t
                        INDEX diff = nki - mki;
                        qgcomplex ival(0,0); 
                        deriv_index_set dset(tindex, t_prime_index, 0, mki, diff, ki, i, a, n, initialguess, regularized);
                        ival += this->trap_rule_integrate(&ssystem::clocal::get_RxPhi_plus_SxG,dset,0,t_prime_index); 
                        // integral from 0 to t'
                        ival += this->trap_rule_integrate(&ssystem::clocal::get_RxPhi,dset,t_prime_index,tindex); 
                        // integral from t' to t
                        sum+=ival; 
                      }
                    }
                  }
                  //Fval = sum; // assign value of integral --> No direct assignments: Negative modes do not get set correctly !!
                  this->assign_reg( Fdest, sum, nki, ki, i, n );
                }}}}
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void calculate_F_G  (unsigned int tindex, unsigned int t_prime_index, bool initialguess, bool nonlinear, regset Fdest)
        {
          if(initialguess){
            Fdest.setprimary();
          }else{
            Fdest.setintermediate();
          }
          // Rs, Ss needed should already be calculated by this point
          INDEX nfs = static_cast<INDEX>(this->numfields);
          int nksize = num_nks; // static_cast<int>( this->nkdomain->size() );
          unsigned int ksize  = num_ks; //static_cast<unsigned int>( this->kdomain_centered->size() );
          // outer indices
          for(int nki = 0; nki < nksize; nki++){
            for(unsigned int ki = 0; ki < ksize; ki++){
              for(INDEX i = 0; i < nfs; i++){
                for(INDEX n = 0; n < nfs; n++){
                  // - - - - - iner contractions (summed over) index loops
                  qgcomplex sum(0,0);
                  for(INDEX a = 0; a < nfs; a++){
                    // do the linear part first ...
                    this->calculate_F_linear_part(tindex, t_prime_index, sum, nki, ki, i, n, a, this->G_data, initialguess);
                    // now do nonlinear part
                    if(nonlinear){
                      for(int mki = -nksize + 1; mki < nksize; mki++){
                        INDEX diff = nki - mki;
                        qgcomplex ival(0,0);
                        deriv_index_set dset( tindex, t_prime_index, t_prime_index, mki, diff, ki, i, a, n, initialguess, false);
                        ival+=trap_rule_integrate(&ssystem::clocal::get_RxG,dset,t_prime_index, tindex);
                        sum+=ival;
                      }
                    }
                  }
                  this->assign_reg(Fdest, sum, nki, ki, i, n );
                  // Fval = sum; // assign value of integral --> No direct assignments, because negative modes do not get written correctly
                }}}} 
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void advance_and_assign(bool initialguess, regset lastTime, 
            regset nextTime, regset Freg, regset Freg_corrector)
        {
          // Note that lastTiime and nextTime will point to the same regset if you are doing an initial guess!!
          const REAL & timestep = this->timestepwidths[this->tindex];
          //const REAL factor = !initialguess ? timestep : 0.5*timestep;
          const REAL factor = 0.5*timestep; // In the case of an initial guess, we march forward 1/2 timestep
          // when not an initial guess, we take the average of two slopes
          INDEX nfs = static_cast<INDEX>(this->numfields);
          unsigned int ksize  = num_ks; //static_cast<unsigned int>( this->kdomain_centered->size() );

          Freg.setprimary(); // Freg is the derivative evaluated at the integer timestep location (the first guess)
          if(!initialguess){
            Freg_corrector.setintermediate(); // Freg_corrected is the derivative evaluated at the half time step (corrector)
          } // else we cannot gaurantee that Freg_corrected != Freg and we would mess it up because they are references!!
          // outer indices
          for(int nki = 0; nki < num_nks; nki++){
            //const int & nk = (*nkdomain)[nki];
            for(unsigned int ki = 0; ki < num_ks; ki++){
              for(INDEX i = 0; i < nfs; i++){
                for(INDEX n = 0; n < nfs; n++){
                  lastTime.setprimary();
                  const qgcomplex last_value = this->get_register(lastTime, nki, ki, i, n);
                  if(initialguess){
                    nextTime.setintermediate(); //  
                    qgcomplex next_value;//       = this->get_register_nonconst(nextTime, nki, ki, i, n);
                    const qgcomplex & F = this->get_register(Freg, nki, ki, i, n);
                    const qgcomplex & vdissip = this->get_register(dissip_register, 0, ki, i, i);
                    REAL dissip = exp( 0.5*timestep*vdissip.real() );
                    next_value = dissip*(last_value + factor*F);
                    this->assign_reg( nextTime, next_value, nki, ki, i, n );
                  }else{
                    nextTime.setprimary(); //
                    qgcomplex next_value;//   = this->get_register_nonconst(nextTime, nki, ki, i, n);
                    const qgcomplex & F1 = this->get_register(Freg, nki, ki, i, n);
                    const qgcomplex & F2 = this->get_register(Freg_corrector, nki, ki, i, n);
                    const qgcomplex & vdissip = this->get_register(dissip_register, 0, ki, i, i);
                    REAL dissip = exp( timestep*vdissip.real() );
                    next_value = dissip*(last_value + factor*(F1+F2)); // factor = 0.5*timestep
                    this->assign_reg( nextTime, next_value, nki, ki, i, n );
                  }
                }}}}
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void advance_and_assign_equal_time(bool initialguess, regset lastTime, regset nextTime, regset Freg_regularized,
            regset Freg_corrector_regularized, regset Freg, regset Freg_corrector)
        {
          const REAL & timestep = this->timestepwidths[this->tindex];
          const REAL factor = 0.5*timestep; // In the case of an initial guess, we march forward 1/2 timestep
          // when not an initial guess, we take the average of two slopes
          INDEX nfs = static_cast<INDEX>(this->numfields);

          Freg.setprimary(); // Freg is the derivative evaluated at the integer timestep location (the first guess)
          if(this->useregularization) Freg_regularized.setprimary();
          if(!initialguess){
            Freg_corrector.setintermediate(); // Freg_corrected is the derivative evaluated at the half time step (corrector)
            if(this->useregularization) Freg_corrector_regularized.setintermediate();
          } // else we cannot gaurantee that Freg_corrected != Freg and we would mess it up because they are references!!
          // outer indices
          for(int nki = 0; nki < num_nks; nki++){
            const int & nk = (nkdomain)[nki];
            for(unsigned int ki = 0; ki < num_ks; ki++){
              for(INDEX i = 0; i < nfs; i++){
                for(INDEX n = 0; n < nfs; n++){
                  lastTime.setprimary();
                  const qgcomplex last_value = this->get_register(lastTime, nki, ki, i, n);
                  if(initialguess){
                    nextTime.setintermediate(); //  
                    qgcomplex next_value;//       = this->get_register_nonconst(nextTime, nki, ki, i, n);
                    qgcomplex F = this->get_register(Freg, nki, ki, i, n);
                    if( this->useregularization ) F += this->get_register(Freg_regularized, nki, ki, i, n);
                    const qgcomplex & vdissip = this->get_register(dissip_register, 0, ki, i, i);
                    REAL dissip = exp( 0.5*timestep*vdissip.real() );
                    qgcomplex Fstar = this->get_register(Freg, nki, ki, n, i);
                    if( this->useregularization ) Fstar += this->get_register(Freg_regularized, nki, ki, n, i);
                    qgcomplex Finhomogeneous = this->get_register( inhomo_forcing_register, nki, ki, i, n );
                    if(nk%2 == 0){
                      next_value = dissip*(last_value + factor*(F+Fstar+Finhomogeneous));
                    }else{
                      next_value = dissip*(last_value + factor*(F-Fstar+Finhomogeneous));
                    }
                    this->assign_reg( nextTime, next_value, nki, ki, i, n );
                  }else{
                    nextTime.setprimary(); //
                    qgcomplex next_value;//   = this->get_register_nonconst(nextTime, nki, ki, i, n);
                    qgcomplex Finhomogeneous = this->get_register( inhomo_forcing_register, nki, ki, i, n );
                    qgcomplex F1 = this->get_register(Freg, nki, ki, i, n);
                    qgcomplex F2 = this->get_register(Freg_corrector, nki, ki, i, n);
                    if( this->useregularization ) {
                      F1 += this->get_register( Freg_regularized, nki, ki, i, n);
                      F2 += this->get_register( Freg_corrector_regularized, nki, ki, i, n);
                    }
                    const qgcomplex & vdissip = this->get_register(dissip_register, 0, ki, i, i);
                    REAL dissip = exp( timestep*vdissip.real() );
                    qgcomplex F1_star = (this->get_register(Freg, nki, ki, n, i));
                    qgcomplex F2_star = (this->get_register(Freg_corrector, nki, ki, n, i));
                    if( this->useregularization ) {
                      F1_star += this->get_register( Freg_regularized, nki, ki, n, i);
                      F2_star += this->get_register( Freg_corrector_regularized, nki, ki, n, i );
                    }
                    if(nk%2 == 0){
                      next_value = dissip*(last_value + factor*( F1 + F1_star + F2 + F2_star + Finhomogeneous)); // factor = 0.5*timestep
                    }else{
                      next_value = dissip*(last_value + factor*( F1 - F1_star + F2 - F2_star +Finhomogeneous )); // factor = 0.5*timestep
                    }
                    this->assign_reg( nextTime, next_value, nki, ki, i, n );
                  }
                }}}}
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void calculate_F_linear_part(unsigned int tindex, unsigned int t_prime_index, qgcomplex & sum, 
            const INDEX nki, const INDEX ki, const INDEX i, const INDEX n, const INDEX a, twotimevar * data, const bool initialguess)
        {
          // Takes Phi or G and returns LijPhi = F_{ij} due to the linear part 
          // this function returns a single value based on the request nki, ki, i, a, n
          // basically return \sum_{npi} L_{ia}^{nki}(ki) * Phi/G_{an}^{nki - npi}(ki)
          // - - - - - iner contractions (summed over) index loops
          // load correct registers
          regset reg = data->getvar(tindex, t_prime_index); 
          if(initialguess) reg.setprimary();
          else reg.setintermediate();
          // sum over angular modes of npi
          const qgcomplex & Ln = this->get_register(lop_register, 0, ki, i, a);
          const qgcomplex & val0 = this->get_register(reg, nki, ki, a, n);
          sum += Ln*val0;
          for(INDEX npi = 1; npi < static_cast<INDEX>(num_nks); npi++){
            const int & np = (nkdomain )[npi];
            const qgcomplex & Lval = this->get_register(lop_register, npi, ki, i, a); 
            const qgcomplex & val1 = this->get_register(reg, nki-npi, ki, a, n);
            const qgcomplex & val2 = this->get_register(reg, nki+npi, ki, a, n);
            qgcomplex lsum = Lval*val1;
            if( np % 2 == 0 ){
              lsum += conj(Lval)*val2;
            }else{
              lsum -= conj(Lval)*val2;
            }
            sum += lsum;
          }
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void copy_regset_to_array(float * data, const int size, const regset & reg, bool real) const
        {
          const int reg_size = num_nks*num_ks*(this->numfields * this->numfields );
          int iter = 0;
          for(INDEX nki = 0; nki < num_nks; nki++){
            for(INDEX ki = 0; ki < num_ks; ki++){
              for(INDEX i = 0; i < static_cast<INDEX>(numfields); i++){
                for(INDEX j = 0; j <static_cast<INDEX>(numfields); j++){
                  rassert( iter < reg_size );
                  if(real){
                    data[iter] = static_cast<float>(this->get_register(reg, nki, ki, i, j).real());
                  }else{
                    data[iter] = static_cast<float>(this->get_register(reg, nki, ki, i, j).imag());
                  }
                  iter++;
                }}}}
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void write_array_to_regset(float * real_data, float * imag_data, const int size, regset reg, bool prov)
        {
          const int reg_size = num_nks*num_ks*(this->numfields * this->numfields );
          if( prov ){
            reg.setintermediate();
          }else{
            reg.setprimary();
          };
          int iter = 0;
          for(INDEX nki = 0; nki < num_nks; nki++){
            for(INDEX ki = 0; ki < num_ks; ki++){
              for(INDEX i = 0; i < static_cast<INDEX>(numfields); i++){
                for(INDEX j = 0; j < static_cast<INDEX>(numfields); j++){
                  rassert( iter < reg_size );
                  this->assign_reg(reg, qgcomplex( real_data[iter], imag_data[iter] ), nki, ki, i, j);
                  iter++;
                }}}}
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        bool save_regset(NcFile * ncfile, const char * name, const twotimevar * reg_data) const
        {
          NcError error_status( NcError::silent_nonfatal );
          NcDim * tdim = ncfile->get_dim("time");
          NcDim * field_dim = ncfile->get_dim("field");
          NcDim * nk_dim = ncfile->get_dim("nkindex");
          NcDim * k_dim = ncfile->get_dim("kindex");

          const NcDim ** dlist = new const NcDim * [6];
          dlist[0] = tdim;
          dlist[1] = tdim;
          dlist[2] = nk_dim;
          dlist[3] = k_dim;
          dlist[4] = field_dim;
          dlist[5] = field_dim;

          string prim_real_name = string(name) + string("_real");
          string prim_imag_name = string(name) + string("_imag");
          string prov_real_name = string(name) + string("_real_provisional");
          string prov_imag_name = string(name) + string("_imag_provisional");

          NcVar * prim_real_var = ncfile->add_var(prim_real_name.c_str(), ncFloat, 6, dlist);
          NcVar * prim_imag_var = ncfile->add_var(prim_imag_name.c_str(), ncFloat, 6, dlist);
          NcVar * prov_real_var = ncfile->add_var(prov_real_name.c_str(), ncFloat, 6, dlist);
          NcVar * prov_imag_var = ncfile->add_var(prov_imag_name.c_str(), ncFloat, 6, dlist);
          if( error_status.get_err() != NC_NOERR ){
            mo.printfmsg("Error writing dissipation operator information.\n");
            mo.printfmsg("NETCDF Error: %s.\n", error_status.get_errmsg());
            return false;
          }
          delete [] dlist;

          prim_real_var->add_att("description", "Real component (t,t',nki,ki,i,j) data.");
          prim_imag_var->add_att("description", "Imaginary component (t,t',nki,ki,i,j) data.");
          prov_real_var->add_att("description", "Real component (t,t',nki,ki,i,j) initial guess data.");
          prov_imag_var->add_att("description", "Imaginary component (t,t',nki,ki,i,j) initial guess data.");

          // int rsize = prim_real_var->num_vals();

          // get the regset size
          const int nk_size = num_nks; //nkdomain->size();
          const int k_size  = num_ks; //kdomain_centered->size();
          const int reg_size = nk_size*k_size*(this->numfields * this->numfields );

          long * edge_lengths = prim_real_var->edges();
          // first check to match
          assert(edge_lengths[2] == nk_size);
          assert(edge_lengths[3] == k_size);
          assert(edge_lengths[4] == static_cast<long>(numfields));
          assert(edge_lengths[5] == static_cast<long>(numfields));

          edge_lengths[0] = 1; // we write 1 time at a time
          edge_lengths[1] = 1; // we write 1 tprime at a time

          long * corner = prim_real_var->edges();
          for(int i = 0; i < 6; i++) corner[i] = 0;

          // allocate space
          float * data = new float[reg_size];

          const unsigned int tsize = this->tindex + 1;
          for(unsigned int ti = 0; ti < tsize; ti++){
            for(unsigned int tpi = 0; tpi <= ti; tpi++){  
              corner[0] = ti;
              corner[1] = tpi;
              prim_real_var->set_cur(corner); prov_real_var->set_cur(corner);
              prim_imag_var->set_cur(corner); prov_imag_var->set_cur(corner);
              //regset cur_reg = (*this.*get_reg)(ti, tpi);
              const regset cur_reg = reg_data->getvar_withcgcmexception(ti,tpi);
              cur_reg.setprimary();
              // copy real primary data first
              copy_regset_to_array( data, reg_size, cur_reg, true );
              if(!prim_real_var->put( data, edge_lengths )){ mo.printfmsg("Failed to write field '%s'.\n", name); return false; };
              float * test = new float[reg_size];
              prim_real_var->get( test, edge_lengths );
              delete [] test;
              copy_regset_to_array( data, reg_size, cur_reg, false );
              if(!prim_imag_var->put( data, edge_lengths )){ mo.printfmsg("Failed to write field '%s'.\n", name); return false; };
              cur_reg.setintermediate();
              copy_regset_to_array( data, reg_size, cur_reg, true );
              if(!prov_real_var->put( data, edge_lengths )){ mo.printfmsg("Failed to write field '%s'.\n", name); return false; };
              copy_regset_to_array( data, reg_size, cur_reg, false );
              if(!prov_imag_var->put( data, edge_lengths )){ mo.printfmsg("Failed to write field '%s'.\n", name); return false; };
            }}
          delete [] edge_lengths;
          delete [] corner;
          delete [] data;	
          return true;
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        bool save_regset_single_time(NcFile * ncfile, const char * name, get_reg_set_single_time_func get_reg) const{
          NcError error_status( NcError::silent_nonfatal );
          NcDim * tdim = ncfile->get_dim("time");
          NcDim * field_dim = ncfile->get_dim("field");
          NcDim * nk_dim = ncfile->get_dim("nkindex");
          NcDim * k_dim = ncfile->get_dim("kindex");

          const NcDim ** dlist = new const NcDim * [5];
          dlist[0] = tdim;
          dlist[1] = nk_dim;
          dlist[2] = k_dim;
          dlist[3] = field_dim;
          dlist[4] = field_dim;

          string prim_real_name = string(name) + string("_real");
          string prim_imag_name = string(name) + string("_imag");
          string prov_real_name = string(name) + string("_real_provisional");
          string prov_imag_name = string(name) + string("_imag_provisional");

          NcVar * prim_real_var = ncfile->add_var(prim_real_name.c_str(), ncFloat, 5, dlist);
          NcVar * prim_imag_var = ncfile->add_var(prim_imag_name.c_str(), ncFloat, 5, dlist);
          NcVar * prov_real_var = ncfile->add_var(prov_real_name.c_str(), ncFloat, 5, dlist);
          NcVar * prov_imag_var = ncfile->add_var(prov_imag_name.c_str(), ncFloat, 5, dlist);
          if( error_status.get_err() != NC_NOERR ){
            mo.printfmsg("Error writing dissipation operator information.\n");
            mo.printfmsg("NETCDF Error: %s.\n", error_status.get_errmsg());
            return false;
          }
          delete [] dlist;

          prim_real_var->add_att("description", "Real component (t,t',nki,ki,i,j) data.");
          prim_imag_var->add_att("description", "Imaginary component (t,t',nki,ki,i,j) data.");
          prov_real_var->add_att("description", "Real component (t,t',nki,ki,i,j) initial guess data.");
          prov_imag_var->add_att("description", "Imaginary component (t,t',nki,ki,i,j) initial guess data.");

          // int rsize = prim_real_var->num_vals();

          // get the regset size
          const int nk_size = num_nks; //nkdomain->size();
          const int k_size  = num_ks; //kdomain_centered->size();
          const int reg_size = nk_size*k_size*(this->numfields * this->numfields );


          long * edge_lengths = prim_real_var->edges();
          // first check to match
          assert(edge_lengths[1] == nk_size);
          assert(edge_lengths[2] == k_size);
          assert(edge_lengths[3] == static_cast<long>(numfields));
          assert(edge_lengths[4] == static_cast<long>(numfields));

          edge_lengths[0] = 1; // we write 1 time at a time

          long * corner = prim_real_var->edges();
          for(int i = 0; i < 5; i++) corner[i] = 0;

          // allocate space
          float * data = new float[reg_size];

          const unsigned int tsize = this->tindex + 1;
          for(unsigned int ti = 0; ti < tsize; ti++){
            corner[0] = ti;
            prim_real_var->set_cur(corner); prov_real_var->set_cur(corner);
            prim_imag_var->set_cur(corner); prov_imag_var->set_cur(corner);
            regset cur_reg = (*this.*get_reg)(ti);
            cur_reg.setprimary();
            // copy real primary data first
            copy_regset_to_array( data, reg_size, cur_reg, true );
            if(!prim_real_var->put( data, edge_lengths )){ mo.printfmsg("Failed to write field '%s'.\n", name); return false; };
            copy_regset_to_array( data, reg_size, cur_reg, false );
            if(!prim_imag_var->put( data, edge_lengths )){ mo.printfmsg("Failed to write field '%s'.\n", name); return false; };
            cur_reg.setintermediate();
            copy_regset_to_array( data, reg_size, cur_reg, true );
            if(!prov_real_var->put( data, edge_lengths )){ mo.printfmsg("Failed to write field '%s'.\n", name); return false; };
            copy_regset_to_array( data, reg_size, cur_reg, false );
            if(!prov_imag_var->put( data, edge_lengths )){ mo.printfmsg("Failed to write field '%s'.\n", name); return false; };
          }
          delete [] edge_lengths;
          delete [] corner;
          delete [] data;	
          return true;
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        bool restore_regset(NcFile * ncfile, const char * name, twotimevar * rdata )
        {
          // restore all time arguments from the file 
          const string vname_real = string(name)  + string("_real");
          const string vname_imag = string(name)  + string("_imag");
          const string vname_real_prov = string(name)  + string("_real_provisional");
          const string vname_imag_prov = string(name)  + string("_imag_provisional");
          NcVar * var_real = ncfile->get_var(vname_real.c_str());
          NcVar * var_imag = ncfile->get_var(vname_imag.c_str());
          NcVar * var_real_prov = ncfile->get_var(vname_real_prov.c_str());
          NcVar * var_imag_prov = ncfile->get_var(vname_imag_prov.c_str());
          if( var_real == NULL || var_imag == NULL || var_real_prov == NULL || var_imag_prov == NULL ){
            fprintf(stderr, "Error: Failed to load two time data set %s from saved model state.\n", name);
            return false;
          }
          const int ndims = var_real -> num_dims();

          assert( ndims == 6 );

          long * edges = var_real->edges();
          edges[0] = 1;
          edges[1] = 1;

          long * corner = var_real->edges();
          for(int i = 0; i < ndims; i++) corner[i] = 0; 

          const int nk_size = num_nks; // nkdomain->size();
          const int k_size  = num_ks; //kdomain_centered->size();
          const int reg_size = nk_size*k_size*(this->numfields * this->numfields );

          float * data_real = new float[reg_size];
          float * data_imag = new float[reg_size];

          for(unsigned int ti = 0; ti <= this->tindex; ti++){
            for(unsigned int tpi = 0; tpi <= ti; tpi++){
              regset reg = rdata->getvar(ti, tpi);
              corner[0] = ti;
              corner[1] = tpi;
              var_real->set_cur(corner); var_real_prov->set_cur(corner);
              var_imag->set_cur(corner); var_imag_prov->set_cur(corner);

              var_real->get( data_real, edges );
              var_imag->get( data_imag, edges );

              write_array_to_regset(data_real, data_imag, reg_size, reg, false); 

              var_real_prov->get(data_real, edges );
              var_imag_prov->get(data_imag, edges );
              write_array_to_regset(data_real, data_imag, reg_size, reg, true); 
            }} 

          delete [] data_real;
          delete [] data_imag;   

          delete [] edges;
          delete [] corner; 
          return true;
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        bool restore_single_time_regset(NcFile * ncfile, const char * name, twotimevar * rdata )
        {
          NcError error_status( NcError::silent_nonfatal );
          // restore all time arguments from the file 
          const string vname_real = string(name)  + string("_real");
          const string vname_imag = string(name)  + string("_imag");
          const string vname_real_prov = string(name)  + string("_real_provisional");
          const string vname_imag_prov = string(name)  + string("_imag_provisional");

          NcVar * var_real = ncfile->get_var(vname_real.c_str());
          NcVar * var_imag = ncfile->get_var(vname_imag.c_str());
          NcVar * var_real_prov = ncfile->get_var(vname_real_prov.c_str());
          NcVar * var_imag_prov = ncfile->get_var(vname_imag_prov.c_str());
          if( var_real == NULL || var_imag == NULL || var_real_prov == NULL || var_imag_prov == NULL ){
            fprintf(stderr, "Error: Failed to load single time data set %s from saved model state.\n", name);
            return false;
          }
          const int ndims = var_real -> num_dims();

          assert( ndims == 5 );

          long * edges = var_real->edges();
          edges[0] = 1;

          long * corner = var_real->edges();
          for(int i = 0; i < ndims; i++) corner[i] = 0; 

          const int nk_size = num_nks; //nkdomain->size();
          const int k_size  = num_ks; //kdomain_centered->size();
          const int reg_size = nk_size*k_size*(this->numfields * this->numfields );

          float * data_real = new float[reg_size];
          float * data_imag = new float[reg_size];

          for(unsigned int ti = 0; ti <= this->tindex; ti++){
            regset reg = rdata->getvar(ti, ti);
            corner[0] = ti;
            var_real->set_cur(corner); var_real_prov->set_cur(corner);
            var_imag->set_cur(corner); var_imag_prov->set_cur(corner);

            var_real->get( data_real, edges );
            var_imag->get( data_imag, edges );

            write_array_to_regset(data_real, data_imag, reg_size, reg, false); 

            var_real_prov->get(data_real, edges );
            var_imag_prov->get(data_imag, edges );
            write_array_to_regset(data_real, data_imag, reg_size, reg, true); 
          } 

          delete [] data_real;
          delete [] data_imag;   

          delete [] edges;
          delete [] corner; 
          return true;
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        void fregdump(FILE * rout, const regset & reginput) const
        {
          regset reg = reginput;
          // uses structured lookup to dump out contents of the register
          size_t ksize = num_ks; //kdomain_centered->size();
          fprintf(rout, "--------------------------------------\n");
          fprintf(rout, "DISPLAYING REGISTER: ");
          fprintf(rout, "@data = %p, @primary = %p, @intermdt = %p\n", reg.get_data(), reg.get_primary_data(),
              reg.get_intermediate_data());
          INDEX nksize = num_nks; //static_cast<INDEX>(nkdomain->size());
          INDEX nfs = static_cast<INDEX>(numfields);
          reg.setprimary();
          fprintf(rout, " ........ %3s %3s %3s %3s\n", "nki", "ki", "i", "n");
          for(INDEX nki = -nksize+1; nki < nksize; ++nki){
            for(INDEX ki = 0; ki < static_cast<INDEX>(ksize); ++ki){
              for(INDEX i = 0; i < nfs; ++i){
                for(INDEX n = 0; n < nfs; ++n){
                  const qgcomplex & val = this->get_register(reg, nki, ki, i, n);
                  fprintf(rout, "  primary %3i %3i %3i %3i = [ %12g, %12g ]\n",
                      static_cast<int>(nki), static_cast<int>(ki), static_cast<int>(i), static_cast<int>(n),
                      val.real(), val.imag());
                }}}}
          reg.setintermediate();
          for(INDEX nki = -nksize+1; nki < nksize; ++nki){
            for(INDEX ki = 0; ki < static_cast<INDEX>(ksize); ++ki){
              for(INDEX i = 0; i < nfs; ++i){
                for(INDEX n = 0; n < nfs; ++n){
                  const qgcomplex & val = this->get_register(reg, nki, ki, i, n);
                  fprintf(rout, " intermdt %3i %3i %3i %3i = [ %12g, %12g ]\n",
                      static_cast<int>(nki), static_cast<int>(ki), static_cast<int>(i), static_cast<int>(n),
                      val.real(), val.imag());
                }}}}
        }
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        }; // end class ssystem::clocal
        /* ------------------------------------------------------------------------------------
           definition of the Ssystem class: horizontally homogeneous PDE solver
           using two-step predictor corrector integration scheme
           -------------------------------------------------------------------------------------- */
        ssystem::ssystem()
        {
          this->local=NULL;
        };
        // ----------------------------------------------------------------------------- //
        ssystem::~ssystem()
        {};
        // ----------------------------------------------------------------------------- //
        void ssystem::setnumbertimesteps(const unsigned int & nsteps){
          local->setnumbertimesteps(nsteps);
        };
        // ----------------------------------------------------------------------------- //
        const REAL & ssystem::gettimestep(const unsigned int & tindexarg) const{
          rassert(tindexarg <= this->local->tindex);
          return this->local->timestepwidths[tindexarg];
        };
        // ----------------------------------------------------------------------------- //
        REAL ssystem::getelapsedtime(const unsigned int & tindexarg) const{
          rassert(tindexarg <= this->local->tindex);
          REAL rvalue = 0;
          for(size_t i = 1; i <= static_cast<size_t>(tindexarg); i++){ // tindex is incremented at the end of a stepforward 
            rvalue += this->local->timestepwidths[i-1];
          }
          return rvalue;
        };
        // ----------------------------------------------------------------------------- //
        void ssystem::sync_tindex(const int ti){
          if( ti >= 0 && ti <= static_cast<int>(this->local->tindex)){
            this->local->tindex = static_cast<unsigned int>(ti);
          }else{
            fatal_stop("Bad argument provided to function.");
          }
        };
        // ----------------------------------------------------------------------------- //
        void ssystem::rewindtime(const unsigned int & nsteps){
          unsigned int & tindex = local->tindex;
          if(nsteps > tindex){
            fatal_stop("Cannot rewind time beyond t = 0.");
          }
          tindex = tindex - nsteps;
        };
        // ----------------------------------------------------------------------------- //
        void ssystem::settprimeintegrationlist(const vector<int> * val){
          this->local->tprimeintegrationlist = new vector<int>;
          *(this->local->tprimeintegrationlist) = *val; // assignment
        }
        // ----------------------------------------------------------------------------- //
        void ssystem::setstepequaltimeonly(bool b){
          this->local->stepequaltimeonly= b;
        };
        // ----------------------------------------------------------------------------- //
        int ssystem::get_tindex() const{
          return this->local->tindex;
        };
        // ----------------------------------------------------------------------------- //
        void ssystem::settimestepwidth(const unsigned int & tindex, const REAL & tstep){
#ifdef DEBUG
          assert(tindex <= this->local->tindex);
          assert(tstep > 0);
#endif
          this->local->timestepwidths[tindex] = tstep;
        };
        // ----------------------------------------------------------------------------- //
        void ssystem::stepforward(REAL timestep, bool initialguess){
          this->local->stepforward(timestep, initialguess);
        }
        // ----------------------------------------------------------------------------- //
        const qgcomplex & ssystem::get_register
          (const regset & reg, const INDEX  m, const INDEX k, const INDEX i, const INDEX n) const
          {
            this->local->get_register_noinline(reg, m, k, i, n);
          }
        // ----------------------------------------------------------------------------- //
        static_container * ssystem::getScoeffs_dyn() {
          return &(local->mo.scoeffs);
        }
        // ----------------------------------------------------------------------------- //
        static_container * ssystem::getRcoeffs_dyn(){
          return &(local->mo.rcoeffs);
        };
        // ------------------------------------------------------------------------------*
        bool ssystem::is_model_using_regularization() const{
          return local->useregularization;
        };
        // ------------------------------------------------------------------------------*
        const regset & ssystem::get_G_set(unsigned int t, unsigned int t_primed) const{
          return this->local->G_data->getvar_withcgcmexception(t,t_primed); 
          // these functions are used by diagnostics
          // blocks so they should throw a cgcmexception if they cannot find what is being requested
        }
        // ------------------------------------------------------------------------------*
        const regset & ssystem::get_F_G_set(unsigned int t, unsigned int tprime) const{
          return this->local->F_G_data->getvar_withcgcmexception(t, tprime);
        };
        // ------------------------------------------------------------------------------*
        const regset & ssystem::get_F_phi_set(unsigned int t, unsigned int tprime) const{
          return this->local->F_phi_data->getvar_withcgcmexception(t, tprime);
        };
        // ------------------------------------------------------------------------------*
        const regset & ssystem::get_S_set(unsigned int t, unsigned int t_primed) const{
          return this->local->S_data->getvar_withcgcmexception(t, t_primed);
        }
        // ------------------------------------------------------------------------------*
        const regset & ssystem::get_R_regularized_set(unsigned int t) const{
          return this->local->R_data_regularized->getvar_withcgcmexception(t,t);
        };
        // ------------------------------------------------------------------------------*
        const regset & ssystem::get_S_regularized_set(unsigned int t) const{
          return this->local->S_data_regularized->getvar_withcgcmexception(t,t);
        };
        // ------------------------------------------------------------------------------*
        const regset & ssystem::get_F_phi_regularized_set(unsigned int t) const{
          return this->local->F_phi_data_regularized->getvar_withcgcmexception(t,t);
        };
        // ------------------------------------------------------------------------------*
        const regset & ssystem::get_R_set(unsigned int t, unsigned int t_primed) const{
          return this->local->R_data->getvar_withcgcmexception(t, t_primed);
        };
        // ----------------------------------------------------------------------------- //
        const regset & ssystem::get_phi_set(unsigned int t, unsigned int t_primed) const{
          return this->local->phi_data->getvar_withcgcmexception(t, t_primed);
        };
	// ----------------------------------------------------------------------------- //
	const twotimevar & ssystem::get_R_data() const {  
		return *(local->phi_data);  
	};
	const twotimevar & ssystem::get_S_data() const {  
		return *(local->S_data);
	};
	const twotimevar & ssystem::get_phi_data() const {  
		return *(local->phi_data); 
	};
	const twotimevar & ssystem::get_G_data() const {   
		return *(local->G_data);
	};
	const twotimevar & ssystem::get_F_phi_data() const {  
		return *(local->F_phi_data);
	};
	const twotimevar & ssystem::get_F_G_data() const {  
		return *(local->G_data);
	};
	const twotimevar & ssystem::get_F_phi_regularized_data() const {  
		return *(local->F_phi_data_regularized);
	};
	const twotimevar & ssystem::get_R_regularized_data() const {  
		return *(local->R_data_regularized);
	};
	const twotimevar & ssystem::get_S_regularized_data() const { 
		return *(local->S_data_regularized);
	};
	// ----------------------------------------------------------------------------- //
	const vector< pair<INDEX, ParserX> > & ssystem::get_dissipation_ops() const {
		return this->local->dissipation_ops; 
	}
	// ----------------------------------------------------------------------------- //
	const vector< pair< rank3 , ParserX> > & ssystem::get_inhomogeneous_ops() const {
		return this->local->inhomogeneous_ops; 
	}
	const lopinfo & ssystem::get_linear_op_info() const {	
		return this->local->linear_ops;
	}
#ifdef DEBUG
	// ----------------------------------------------------------------------------- //
	void ssystem::dumpFPhi_regularized(unsigned int t) const{
		FILE * rout = getrout();
		fprintf(rout, "Retrieve FPhi regularized data for t=%u\n", t);
		try {
			const regset & reg = this->get_F_phi_regularized_set(t);
			this->regdump(reg);
		}catch(cgcmexception & e){
			fprintf(rout, "%s\n", e.what());
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dumpR_regularized(unsigned int t) const{
		FILE * rout = getrout();
		fprintf(rout, "Retrieve R regularized data for t=%u\n", t);
		try {
			const regset & reg = this->get_R_regularized_set(t);
			this->regdump(reg);
		}catch(cgcmexception & e){
			fprintf(rout, "%s\n", e.what());
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dumpS_regularized(unsigned int t) const{
		FILE * rout = getrout();
		fprintf(rout, "Retrieve S regularized data for t=%u\n", t);
		try {
			const regset & reg = this->get_S_regularized_set(t);
			this->regdump(reg);
		}catch(cgcmexception & e){
			fprintf(rout, "%s\n", e.what());
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dumpR(unsigned int t, unsigned int t2) const{
		FILE * rout = getrout();
		fprintf(rout, "Retrieve R data for t1=%u, t2=%u\n", t,t2);
		try { 
			const regset & reg = this->get_R_set(t,t2);
			this->regdump(reg);
		}catch(cgcmexception & e){
			fprintf(rout, "%s\n", e.what());
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dumpS(unsigned int t, unsigned int t2) const{
		FILE * rout = getrout();
		fprintf(rout, "Retrieve S data for t1=%u, t2=%u\n", t,t2);
		try { 
			const regset & reg = this->get_S_set(t,t2);
			this->regdump(reg);
		}catch(cgcmexception & e){
			fprintf(rout, "%s\n", e.what());
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dumpPhi(unsigned int t, unsigned int t2) const{
		FILE * rout = getrout();
		fprintf(rout, "Retrieve Phi data for t1=%u, t2=%u\n", t,t2);
		try {
			const regset & reg = this->get_phi_set(t,t2);
			this->regdump(reg);
		}catch(cgcmexception & e){
			fprintf(rout, "%s\n", e.what());
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dumpG(unsigned int t, unsigned int t2) const{
		FILE * rout = getrout();
		fprintf(rout, "Retrieve G data for t1=%u, t2=%u\n", t,t2);
		try {
			const regset & reg = this->get_G_set(t,t2);
			this->regdump(reg);
		}catch(cgcmexception & e){
			fprintf(rout, "%s\n", e.what());
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dumpFG(unsigned int t, unsigned int t2) const{
		FILE * rout = getrout();
		fprintf(rout, "Retrieve FG data for t1=%u, t2=%u\n", t,t2);
		try {
			const regset & reg = this->get_F_G_set(t,t2);
			this->regdump(reg);
		}catch(cgcmexception & e){
			fprintf(rout, "%s\n", e.what());
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dumpFPhi(unsigned int t, unsigned int t2) const{
		FILE * rout = getrout();
		fprintf(rout, "Retrieve FPhi data for t1=%u, t2=%u\n", t,t2);
		try {
			const regset & reg = this->get_F_phi_set(t,t2);
			this->regdump(reg);
		}catch(cgcmexception & e){
			fprintf(rout, "%s\n", e.what());
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dump_all_FPhi_regularized() const{
		for(unsigned int ti = 0; ti <= local->tindex; ti++){
			this->dumpFPhi_regularized(ti);
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dump_all_R_regularized() const{
		for(unsigned int ti = 0; ti <= local->tindex; ti++){
			this->dumpR_regularized(ti);
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dump_all_S_regularized() const{
		for(unsigned int ti = 0; ti <= local->tindex; ti++){
			this->dumpS_regularized(ti);
		}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dump_all_S() const{
		for(unsigned int ti = 0; ti <= local->tindex; ti++){
			for(unsigned int tpi = 0; tpi <= ti; tpi++){
				this->dumpS(ti, tpi);
			}}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dump_all_R() const{
		for(unsigned int ti = 0; ti <= local->tindex; ti++){
			for(unsigned int tpi = 0; tpi <= ti; tpi++){
				this->dumpR(ti, tpi);
			}}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dump_all_G() const{
		for(unsigned int ti = 0; ti <= local->tindex; ti++){
			for(unsigned int tpi = 0; tpi <= ti; tpi++){
				this->dumpG(ti, tpi);
			}}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dump_all_FG() const{
		for(unsigned int ti = 0; ti <= local->tindex; ti++){
			for(unsigned int tpi = 0; tpi <= ti; tpi++){
				this->dumpFG(ti, tpi);
			}}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dump_all_FPhi() const{
		for(unsigned int ti = 0; ti <= local->tindex; ti++){
			for(unsigned int tpi = 0; tpi <= ti; tpi++){
				this->dumpFPhi(ti, tpi);
			}}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::dump_all_Phi() const{
		for(unsigned int ti = 0; ti <= local->tindex; ti++){
			for(unsigned int tpi = 0; tpi <= ti; tpi++){
				this->dumpPhi(ti, tpi);
			}}
	};
	// ----------------------------------------------------------------------------- //
	void ssystem::regdump(const regset & reginput) const{
		this->fregdump(getrout(), reginput);
	};
#endif
	// ----------------------------------------------------------------------------- //
	void ssystem::fregdump(FILE * rout, const regset & reginput) const{
		regset reg = reginput;
		// uses structured lookup to dump out contents of the register
		size_t ksize = local->num_ks; //kdomain_centered->size();
		fprintf(rout, "--------------------------------------\n");
		fprintf(rout, "DISPLAYING REGISTER: ");
		fprintf(rout, "@data = %p, @primary = %p, @intermdt = %p\n", reg.get_data(), reg.get_primary_data(),
				reg.get_intermediate_data());
		INDEX nksize = static_cast<INDEX>(local->nkdomain.size());
		INDEX nfs = static_cast<INDEX>(local->numfields);
		reg.setprimary();
		fprintf(rout, " ........ %3s %3s %3s %3s\n", "nki", "ki", "i", "n");
		for(INDEX nki = -nksize+1; nki < nksize; ++nki){
			for(INDEX ki = 0; ki < static_cast<INDEX>(ksize); ++ki){
				for(INDEX i = 0; i < nfs; ++i){
					for(INDEX n = 0; n < nfs; ++n){
						const qgcomplex & val = this->get_register(reg, nki, ki, i, n);
						fprintf(rout, "  primary %3i %3i %3i %3i = [ %12g, %12g ]\n",
								static_cast<int>(nki), static_cast<int>(ki), static_cast<int>(i), static_cast<int>(n),
								val.real(), val.imag());
					}}}}
		reg.setintermediate();
		for(INDEX nki = -nksize+1; nki < nksize; ++nki){
			for(INDEX ki = 0; ki < static_cast<INDEX>(ksize); ++ki){
				for(INDEX i = 0; i < nfs; ++i){
					for(INDEX n = 0; n < nfs; ++n){
						const qgcomplex & val = this->get_register(reg, nki, ki, i, n);
						fprintf(rout, " intermdt %3i %3i %3i %3i = [ %12g, %12g ]\n",
								static_cast<int>(nki), static_cast<int>(ki), static_cast<int>(i), static_cast<int>(n),
								val.real(), val.imag());
					}}}}
	};
	// ----------------------------------------------------------------------------- //
	bool ssystem::restore_model_state(NcFile * ncfile){
		// load the dissipation operators
		try {
			check_dissip_ops( ncfile, this->local->dissipation_ops );
			check_inhomogeneous_ops( ncfile, this->local->inhomogeneous_ops );
			check_linear_ops ( ncfile, this->local->linear_ops.lop_list );

			bool nonlinear = local->mo.get_global_param_by_name(gdb_str_NONLINEAR).getiValue() != 0;

			NcVar * vtimesteps = ncfile->get_var("timesteps");
			if( vtimesteps == NULL ){ // we have a problem
				fprintf(stderr, "Error: Could not find variable '%s' in model saved state.\n", "timesteps");
				fprintf(stderr, "       Is model saved state file corrupt?\n");
				return false;
			}
			// adjust the time window
			NcDim * tdim = ncfile->get_dim("time"); // dimension has been retrieved before, so if it failed it should have halted
			const int ntimes = static_cast<int>(tdim -> size()); // the running of the model already
			float * tsteps = new float[ntimes];
			if(!vtimesteps->get(tsteps, ntimes, 0, 0, 0, 0)){
				fprintf(stderr, "Error: Could not load variable '%s' in model saved state.\n", "timesteps");
				fprintf(stderr, "       Is model saved state file corrupt?\n");
				return false;
			}
			if( this->local->tindex < static_cast<unsigned int>(ntimes - 1) ){
#ifndef DEBUG
				delete [] this->local->timestepwidths;
				local->timestepwidths = NULL;
#endif
				this->setnumbertimesteps( ntimes );
				// now we need to add time if needed
				// fill in previous ones
				while(local->tindex < static_cast<unsigned int>(ntimes) ){
					local->addtime(local->tindex, tsteps[local->tindex], true); // we need a true initialguess to cause 
					local->tindex++;  							  // an actual expansion of the two time reg set
				};
				local->tindex = ntimes - 1; 
			}

			local->restore_regset( ncfile, "phi", local->phi_data); 
			local->restore_regset( ncfile, "Fphi", local->F_phi_data); 
			local->restore_regset( ncfile, "G", local->G_data); 
			local->restore_regset( ncfile, "FG", local->F_G_data); 
			local->restore_regset( ncfile, "R", local->R_data); 
			local->restore_regset( ncfile, "S", local->S_data);

			if( local->useregularization ){
				// make sure that regularized registers exists
				if(!local->restore_regset( ncfile, "Fphi_regularized", local->F_phi_data_regularized)) return false;
				if(!local->restore_regset( ncfile, "R_regularized", local->R_data_regularized)) return false;
				if(!local->restore_regset( ncfile, "S_regularized", local->S_data_regularized)) return false;
			}

			// now we need to restore the registers .....
			return true;
		}catch( bool ){ 
			return false;
		}
	};
	// ----------------------------------------------------------------------------- //
	bool ssystem::save_model_state(NcFile * ncfile) const{
		// save dissipation operators first
          model_class & mo = local->mo;
          NcError error_status( NcError::silent_nonfatal );
          const vector< pair<INDEX, ParserX> > & dissipation_ops = local->dissipation_ops;

          const int num_d_ops = static_cast<int>(dissipation_ops.size());
          for(int iter = 0; iter < num_d_ops; iter++){
            string name("dissip_op_");
            name+=convertint2string(iter);
            NcVar * var = ncfile->add_var(name.c_str(), ncInt, NULL);
            var->add_att("module", "ssystem");
            var->add_att("type", "dissipation operator");
            const int field_index = static_cast<int>(dissipation_ops[iter].first);
            var->add_att("field_index", static_cast<int>(field_index));
            var->add_att("expr", dissipation_ops[iter].second.GetExprAsString().c_str());
            var->put( & field_index, 1, 0,0,0,0 );
            if( error_status.get_err() != NC_NOERR ){
              mo.printfmsg("Error writing dissipation operator information.\n");
              mo.printfmsg("NETCDF Error: %s.\n", error_status.get_errmsg());
              return false;
            }
          }
          const vector< pair<rank3, ParserX> > & inhomogeneous_ops = local->inhomogeneous_ops;
          // save inhomogeneous operators
          const int num_inh_ops = static_cast<int>(inhomogeneous_ops.size());
          for(int iter = 0; iter < num_inh_ops; iter++){
            string name("inhomogeneous_op_");
            name+=convertint2string(iter);
            NcVar * var = ncfile->add_var(name.c_str(), ncInt, NULL);
            var->add_att("module", "ssystem");
            var->add_att("type", "inhomogeneous operator");
            const rank3 & r3 =  inhomogeneous_ops[iter].first;
            var->add_att("nkindex", static_cast<int>(r3[0]) );
            var->add_att("Iindex", static_cast<int>(r3[1]) );
            var->add_att("Jindex", static_cast<int>(r3[2]) );
            var->add_att("expr", inhomogeneous_ops[iter].second.GetExprAsString().c_str());
            if( error_status.get_err() != NC_NOERR ){
              mo.printfmsg("Error writing inhomogeneous operator information.\n");
              mo.printfmsg("NETCDF Error: %s.\n", error_status.get_errmsg());
              return false;
            }
          }
          const lopinfo & linear_ops = local->linear_ops;
          // save linear operator information (fun fun)
          const int num_lops = static_cast<int>(linear_ops.lop_list.size());
          for(int iter = 0; iter < num_lops; iter++){
            string name("linear_op_"); name+=convertint2string(iter);
            const lopinfo::element & lop = linear_ops.lop_list[iter];
            NcVar * var = ncfile->add_var(name.c_str(), ncInt, NULL);
            var->add_att("module", "psystem");
            var->add_att("type", "linear operator");
            var->add_att("nkindex", static_cast<int>(lop.angular_mode) ); 
            var->add_att("Iindex", static_cast<int>(lop.field_indices[0]));
            var->add_att("Jindex", static_cast<int>(lop.field_indices[1]));
            var->add_att("expr", lop.parser->GetExprAsString().c_str());
            if( error_status.get_err() != NC_NOERR ){
              mo.printfmsg("Error writing linear operator information.\n");
              mo.printfmsg("NETCDF Error: %s.\n", error_status.get_errmsg());
              return false;
            }
          }
          // save the TIMESTEP size data
          const unsigned int tsize = this->local->tindex+1;
          if(tsize > 0){
            NcDim * tdim = ncfile->get_dim("time");
            { NcVar * var = ncfile->add_var("timesteps", ncFloat, tdim);
              float * data = new float[tsize]; 
              for(unsigned int iter = 0; iter < tsize; iter++) data[iter] = static_cast<float>(this->local->timestepwidths[iter]);
              data[local->tindex] = static_cast<float>(mo.get_global_param_by_name(gdb_str_TIMESTEP).getfValue());
              var->add_att("description", "The nondimensional timestep at each time step.");
              var->put(data, tsize, 0, 0, 0, 0);
              delete [] data;
            }

            if( ! local->save_regset(ncfile, "phi", local-> phi_data) ) return false; 
            if( ! local->save_regset(ncfile, "G", local-> G_data) ) return false; 
            if( ! local->save_regset(ncfile, "Fphi", local-> F_phi_data) ) return false; 
            if( ! local->save_regset(ncfile, "FG", local-> F_G_data) ) return false; 
            if( ! local->save_regset(ncfile, "R", local-> R_data) ) return false; 
            if( ! local->save_regset(ncfile, "S", local-> S_data) ) return false; 
            if(local->useregularization){
              if( ! local->save_regset(ncfile, "Fphi_regularized", local-> F_phi_data_regularized) ) return false; 
              if( ! local->save_regset(ncfile, "R_regularized", local-> R_data_regularized) ) return false; 
              if( ! local->save_regset(ncfile, "S_regularized", local-> S_data_regularized) ) return false; 
            }
          } 
          return true;
        };
        // ----------------------------------------------------------------------------- //
        REAL ssystem::getflux(unsigned int i, unsigned int j, unsigned int tindex, bool initial_guess) const
        {
          return local->getflux(i,j,tindex, initial_guess);
        }
        // ----------------------------------------------------------------------------- //
        void ssystem::dumpdata() const{ 
          local->dumpdata();

        }
        // ----------------------------------------------------------------------------- //
        void ssystem::init(model_class & mo, spreprocdata & spdata){

          bool nonlinear = mo.get_global_param_by_name(gdb_str_NONLINEAR).getiValue() != 0;

          bool useregularization = false; // this is set to true if slot regularization is supposed to occur
          if(nonlinear && (mo.get_global_param_by_name(gdb_str_IGNORETRIANGLEBASE).getiValue() != 0)){
            useregularization = true;
          }

          // create the local class
          this->local = new clocal( mo, *this, static_cast<int>(mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue()),
              mo.get_global_prof_by_name(gdb_prof_str_NK_DOMAIN).getiData(),
              mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData(),
              mo.get_global_prof_by_name(gdb_prof_str_DK_DOMAIN).getData(),
              mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED_SQUARED).getData(),
              nonlinear,
              useregularization
              );
          if( ! this->local->init(spdata) ){
            throw cgcmexception("S-system failed to initialize.", -1);
          }
        }


